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ABSTRACT 

Wc discuss a hydrodynamical model for the dispersal of protoplanetary discs around 
young, low mass (< 1.5 Mq) stars by photoevaporation from the central object's 
energetic radiation, which considers the far-ultraviolet as well as the X-ray component 
of the radiation field. We present analytical scaling relations and derive estimates for 
the total mass-loss rates, as well as discussing the existence of similarity solutions 
for flows from primordial discs and discs with inner holes. Furthermore, we perform 
numerical calculations, which span a wide range of parameter space and allow us to 
provide accurate scalings of the mass-loss rates with the physical parameters of the 
systems (X-ray and FUV luminosity, stellar mass, disc mass, disc temperature and 
inner hole radius). 

The model suggest that the X-ray component dominates the photoevaporative 
mass-loss rates from the inner disc. The mass-loss rates have values in the range from 
10""'^"'^ to 10~^ Mq yr~-'^and scale linearly with X-ray luminosity, with only a weak 
dependence on the other parameters explored. However, in the case of high FUV to 
X-ray {Lpuy / Lx > 100) luminosity ratios, the FUV constricts the X-ray flow and 
may dominate the mass-loss. 

Simulations of low mass discs with inner holes demonstrate a further stage of disc 
clearing, which we call 'thermal sweeping'. This process occurs when the mid-plane 
pressure drops to sufficiently low values. At this stage a bound, warm. X-ray heated 
region becomes sufficiently large and unstable, such that the remaining disc material 
is cleared on approximately dynamical time-scales. This process significantly reduces 
the time taken to clear the outer regions of the disc, resulting in an expected transition 
disc population that will be dominated by accreting objects, as indicated by recent 
observations. 
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1 INTRODUCTION 

The evolution and dispersal of protoplanetary discs is an 
important step in the star and planet formation process. 
Protoplanetary discs provide the material from which plan- 
ets form, implying that the initial evolution of newly formed 
planets is coupled to the properties of the parent disc. More- 
over the disc dispersal time-scale ultimately sets the time in 
which (gas) planets must form. 

It is observationally well established that at an age of 
1 Myr most stars are surrounded by discs that are optically 
thick at infra-red (IR) wavelengths, which have in most cases 
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disappeared by typical ages of 10 Myr (e.g. Haisch et al. 
2001, Mamajek 2009). Optically thick discs around young 
stars evolve in a manner that is consistent with 'standard' 
viscous evolution theory, where the disc diagnostics decline 
with time in an approximately power law fashion (Hartmann 
et al. 1998). However, the final evolution from disc-bearing 
(primordial) to disc-less status appears to be much more 
rapid than viscous theory predicts (e.g Luhman et al. 2010; 
MuzeroUe et al. 2010; Ercolano et al. 2011) and to occur 
over the entire range of disc radii probed by IR observa- 
tions. Furthermore, given that the sub-millimetre observa- 
tions show that most non-accreting stars (WTTs) are de- 
void of emission out to several hundreds of AU (Duvert et 
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al. 2000, Andrews & Williams 2005, 2007), the clearing must 
be correlated across all disc radii from lAU to lOOAU. 

Further clues regarding the mechanism for clearing discs 
are provided by objects which have been identified as 'tran- 
sition' discs. These are discs which show a deficit in opac- 
ity at near-IR (NIR) wavelengths compared to a standard 
optically thick 'primordial' disc, but are consistent with 
standard optically thick discs at Mid-IR/Far-IR (MIR/FIR) 
wavelengths; furthermore some of these objects show small 
NIR excesses above the photosphere (not consistent with a 
primordial disc) have been labelled as 'pre-transition' discs 
(e.g. Espaillat et al. 2010; Furlan et al. 2011). Such spectral 
characteristics are most readily explained in terms of a hole 
or gap in the dust in the inner regions of the disc. However, 
there is no clear consensus as to the origin of such structures 
and the debate is complicated by the fact that they display a 
wide variety of other observational characteristics: for exam- 
ple some show evidence for accretion (e.g. Calvet et al. 2005, 
Najita et al. 2007, Hughes et al. 2009) while others do not 
(Cieza et al. 2010; Merin et al. 2010); there is also evidence 
in some systems that the deficit in opacity at NIR wave- 
lengths does not necessarily preclude significant quantities 
of gas and small amounts of dust in the inner disc (Espaillat 
et al. 2010) although this is not always the case (Calvet et 
al. 2002). The frequency of these 'transition/pre-transition' 
discs compared to those of primordial discs suggest that ap- 
proximately 10-20% of discs are in this stage (Strom et al. 
1989, Skrutskie et al. 1990, Luhman et al. 2010, Ercolano 
et al. 2011a; Furlan et al. 2011). Assuming that these ob- 
jects are in transition from a primordial to a fully cleared 
state, this observation provides two important constraints: 
the transition time from primordial to cleared is approxi- 
mately 10% of the discs total lifetime, and that the clearing 
process proceeds from the inside out. 

While many theories have been developed in order to 
explain the observations of 'transition' discs: grain growth 
(DuUemond & Dominiki, 2005); photophoresis (Krasuss et 
al. 2007); MRI-driven winds (Suzuki & Inutsuka, 2009), only 
two mechanism have been proposed that can explain the ob- 
servations of 'transition' discs in the context of total (gas 
& dust) disc clearing (assuming that indeed all discs pass 
through a 'transition' disc phase). These are photoevapo- 
ration (e.g. Clarke et al. 2001) and planet formation (e.g. 
Armitage & Hansen, 1999). Photoevaporation and planet 
formation are processes which compete to remove gas from 
the disc, and it is still a matter for debate whether one pro- 
cess dominates over the other. Ercolano & Clarke (2010) 
discuss the implications of metallicity on both processes and 
Yasui et al. (2009) provide very preliminary observations of 
smaller disc frequencies in low metallicity regions in the ex- 
treme outer Galaxy, consistent with the predictions from 
X-ray photoevaporation models and arguing against planet 
formation as a main dispersal mechanism. 

Like any theory, the photoevaporation model has pro- 
gressed through several stages of evolution, with the orig- 
inal theories only including photoevaporation without vis- 
cous disc evolution (Shu et al. 1993; HoUenbach et al. 1994; 
Yorke & Welz 1996; Richling & Yorke 1997), or photoevapo- 
ration from external sources (Johnstone et al. 1998; Richling 
& Yorke 1998; Richling & Yorke 2000; Adams et al. 2004). It 
was Clarke et al. (2001) who realised that mass-loss from the 
surface of a disc due to heating (photoevaporation) , coupled 



with viscous evolution results in a gap opening in the in- 
ner disc when the photoevaporation rate and accretion rate 
become comparable. The inner disc then drains onto the 
central star on its own much faster local viscous time-scale. 
While the basic ideas have remained identical (photoevapo- 
ration and viscous process compete until photoevaporation 
eventually opens a gap in the inner disc which then rapidly 
drains on to the central star), the actual calculations of the 
photoevaporation rates have progressed greatly. The origi- 
nal model used semi-analytic mass-loss rates due to EUV 
radiation calculated by HoUenbach et al. (1994), while later 
the photoevaporative fiow was modelled hydrodynamically 
by Font et al. (2004). In these models the disc is mainly 
irradiated by a diffuse field of recombination photons from 
the atmosphere of the inner disc. Subsequently Alexander et 
al. (2006a) realised that, following the draining of the inner 
disc, the inner-edge of the hole is instead directly exposed 
to the radiation field from the star; this results in higher 
mass-loss rates which rapidly clear the remaining disc out 
to large radii (Alexander et al. 2006b). 

Recently, a lot of attention has been paid to other high 
energy radiation fields; X-rays (Alexander et al. 2004, Er- 
colano et al. 2008,2009, Owen et al. 2010,2011b) and FUV 
(Gorti & HoUenbach 2008,2009; Gorti et al. 2009). These cal- 
culations have yielded much higher photoevaporation rates 
than the original EUV rate. These higher mass-loss rates pri- 
marily arise from significant mass-loss at larger radii in the 
disc {R >10 AU), due to the ability of the FUV and X-rays 
to heat much larger columns of gas to of order the escape 
temperature. While the central star certainly produce signif- 
icant fiuxes of FUV, EUV and X-rays (e.g. Alexander et al. 
2005, Ingleby et al. 2011) it is important to determine which, 
or which combination of high energy rates are responsible 
for setting the photoevaporation rate. 

Some of the main observational tests of photoevapora- 
tion are based on observational statistics of disc populations 
as a function of time. Owen et al. (2011b) showed that the 
X-ray photoevaporation model is consistent with the cur- 
rent diagnostics of disc evolution, and made some predic- 
tions about how these diagnostics should correlate with X- 
ray luminosity. Furthermore, Drake et al. (2009) reported 
that accretion rate measurements and X-ray observations in 
Orion, point to a process of X-ray driven 'photoevaporation- 
starved accretion', where higher mass- loss rates prevent ma- 
terial accreting onto the central star, leading to lower accre- 
tion rates. This process was theoretically demonstrated by 
Owen et al. (2011b) and is also likely to explain the system- 
atic discrepancy between the X-ray luminosities of accreting 
and non-accreting young stars. Unlike the previous EUV 
models (e.g. Alexander et al. 2006b; Alexander & Armitage, 
2009), the X-ray photoevaporation model can account for 
a large fraction of the observed 'transition' discs. In par- 
ticular, photoevaporation can account for the population of 
small inner holes R < 20 AU, those with low accretion rates 
Mt < 1O~*M0 yr~^and those with gas inside the inner-hole. 
However, 'transition' discs with large inner holes and large 
accretion rates are likely to be created through a different 
mechanism (e.g. Espaillat et al. 2010). 

The first direct evidence and best test of individual 
photoevaporation models is the detection of a blue-shifted 
12.8/im Nell line from the surface of the disc around TW 
Hya (Pascucci & Sterzik 2009), which is consistent with ei- 
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ther an EUV or X-ray driven photoevaporative wind (Pas- 
cucci et al. 2011). Moreover, Ercolano & Owen (2010) re- 
cently showed that the previously unexplained low veloc- 
ity component of the OI 6300A, - emission blue-shifted by 
~ 5km detected by Hartigan et al. 1995 around most 
primordial discs - is naturally produced in an X-ray driven 
photoevaporative flow. The observed luminosity is too high 
to be consistent with an EUV wind and the blue-shift is 
too large to be consistent with an FUV wind. Currently, 
the only well resolved OI 6300A line is detected from the 
accreting 'transition' disc around TW Hya, where the line 
presents with no blue-shift (Pascucci et al. 2011). Unfortu- 
nately, hydrodynamic models of the inner regions of photoe- 
vaporating transition discs during the accreting phase do not 
exist for any of the photoevaporation models (only transi- 
tion discs with cleared inner gas discs were considered by Er- 
colano & Owen 2010) and useful comparisons must wait until 
these calculations are completed. Therefore, pending these 
calculations, and the availability of high resolution spectro- 
scopic observations of many objects, direct comparisons of 
predicted to observed OI line profiles are not possible. 

While observations are providing a fruitful avenue for 
direct comparisons of photoevaporation models, in this work 
we develop a theoretical framework in which we can anal- 
yse the combined effects off all three heating mechanisms 
(EUV, FUV and X-rays). The most relevant scale in any pho- 
toevaporation model is the radius at which a gas of a given 
temperature (T^as) becomes unbound from the central star 
and is given bjo: 

'^.-8.9AU(^)'^(J^) (1) 

For the EUV photoevaporation model this becomes a fixed 
radius (for a given disc), since EUV heating gives rise to 
essentially isothermal gas at 10*K. For the X-ray and FUV 
model it is better to think of this radius as an escape tem- 
perature scaled length, allowing the consideration of differ- 
ent radial scales around different mass stars. However, in 
order to be consistent with the previous EUV only litera- 
ture (HoUenbach et al. 1994; Clarke et al. 2001; Font et 
al. 2004; Alexander et al. 2006a), we separately fix Rg = 
8.9(1/1M0)AU. Furthermore, as we will discuss in detail in 
Section 2, this radius has an important physical significance 
as it is a measure of the radius at which a flow at tempera- 
ture Tgas can go through the sonic surface. For X-ray heating 
in the range Tgas ~ 1000- lOOOOK, rg ~ 1- lOAU for a very 
low mass star (O.IM©) and rg ~ 10— lOOAU for a solar type 
star. For FUV heating, Tgas < lOOOK, > lOAU for a very 
low mass star and rg > lOOAU for solar type stars. 

In this work, we lay out the theoretical foundation for 
photoevaporation, deriving the scaling with X-ray luminos- 
ity and stellar mass in Sections 2. In Section 3 we discuss 
the likely interplay between the effects of combined EUV, 
FUV and X-ray heating on the calculation of photoevap- 
oration rates. Section 4 presents the results of radiation- 
hydrodynamic calculations that test the theoretical founda- 
tions laid out in the Section 2, while Section 5 contains some 

^ Throughout this work we use {R, 4>, z} to represent the cylin- 
drical co-ordinate system and {r,8,(p} to represent the spherical 
co-ordinate system, in all cases centred on the star. 



simple numerical experiments that address the discussion of 
combined FUV and X-ray heating contained in Section 3. 
Section 6 describes the final clearing of the outer disc by 
'thermal sweeping'. In Section 7 we describe the implica- 
tions of photoevaporation on disc evolution and dispersal. 
In Section 8 we discuss the presented results and draw our 
conclusions in Section 9. 

In Appendix A we discuss the possibility of self- 
similarity in X-ray heated flows and in Appendix B we pro- 
vide numerical flts to our photoevaporation calculations for 
use by the community. 



2 THE THEORY OF X-RAY 
PHOTOEVAPORATION 

The problem of thermally driven winds from discs has been 
studied analytically in the past (Begelman et al. 1983; Liff- 
man 2003; Adams et al. 2004). While a full analytic solution 
to the problem has not been obtained in even the most sim- 
plistic thermodynamic constructions, considerable progress 
can still be made. In the Appendix we discuss the existence 
of similarity solutions with regard to the problem of X-ray 
photoevaporation. The results from the analysis are as fol- 
lows: 

• For primordial discs, the total mass-loss rates scale lin- 
early with X-ray luminosity and are independent of stellar 
mass. 

• For primordial discs, the density at a given value of 
r/Rg scales with X-ray luminosity and stellar mass as 
n{Lx, Mr) oc Lx M^"^ for primordial discs. 

• For discs with inner holes of fixed radius, the total mass- 
loss scales linearly with the X-ray luminosity. 

• Discs with inner holes with a range of radii are not 
strictly mutually self-similar: see section 4.3.1 for an ap- 
proximate self-similarity argument that can be used to un- 
derstand the numerical results in this case. 

In this section, we derive the scalings listed above using 
a series of simple analytic arguments. Since the mass-flux 
is conserved in any steady-state photoevaporative flow we 
are free to evaluate this quantity at any point in the flow. 
Therefore, choosing to do this at the sonic surface the total 
mass-loss rate is given by: 

M„ = / n,(rs)c,(r,)l(rs)-dS (2) 
J s 

where S is the sonic surface, & Ca are the density and 
sound speed at the sonic surface respectfully, and 1 is the 
unit vector along a streamline. In the case of a flow from a 
disc, this can be approximately written as: 

/■oo 

M» = 2/ 2'KRns{R)cs{R)AR (3) 
Jo 

where we note that the true mass loss flux will differ from 
this by a factor of order unity due to the fact that the flow is 
not purely vertical at the sonic surface. Therefore, provided 
we can determine the flow properties at the sonic surface we 
can estimate the total mass-loss rates. 
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2.1 The Sonic Surface in X-ray Heated Flows 

The sonic point is a critical point of the steady-state mo- 
mentum equation when u = Cs, and physically arises as a 
consequence of wave-steeping requirements, namely that no 
strong shock forms (Landau & Lifshitz, 1987). 

Several authors have derived the conditions on the 
sound speed at the sonic surface for a flow emerging from 
a disc (e.g. Begelman et al. 1983; Adams et al. 2004; also 
shown in Appendix A for completeness). Under the con- 
ditions expected for disc winds (sub-sonic launching, flow 
time-scale > thermal time-scale, height of the sonic surface 
[zs) ^ R) the sound speed at the sonic surface should be 
to order unity the Parker value at each cylindrical radius, 
namely: 



2R 



(4) 



Therefore, we can a priori determine the sound speed, and 
hence the temperature at the sonic surface. 

2.1.1 Density- Sound- Speed Coupling in X-ray Heated 
Flows 

If we can specify both the density and sound speed at the 
sonic surface then we have determined the mass-flux in the 
entire flow as mass-flux is conserved. We argued above that 
in the case of disc winds, the sound speed at the sonic sur- 
face is roughly given by the Parker value (Equation |4]). This 
means that the gas temperature at the sonic surface is also 
simply given as a function of cylindrical radii. 

It is well known that the temperature of X-ray heated 
gas can be described in terms of the local ionization param- 
eter (e.g. Tarter et al. 1969), given by: 



AttFx 



(5) 



where Fx is the X-ray flux (Lx/47rr'^). Ercolano et al. 
(2009ab) & Owen et al. (2010) found that the temperature 
of gas within a column density to the central star of lO^'^ 
cm~^ (this being the penetration depth of 1 KeV photon^ 
to the central star can be described in terms of optically 
thin X-ray heating. Thus, adopting the assumption of opti- 
cally thin X-ray heating the ionization parameter becomes 
f = Lxjnr^ , and the gas temperature is specifled in terms of 
local variables only, where the form of this profile T — /(^) 
is shown in Figure [1] 

At the sonic surface (provided ^ i? so i? ~ r) we 
may write: 



fcs/ {Lx/usR^) GAh 



2R 



(6) 



Therefore, we may re-arrange Equationj6]to give an expres- 
sion for the density at the sonic surfac4f|: 



Lx_ 



1 / GM,^mh \ 

^ \ 2kBR ) 



(7) 



^ see discussion in Section 4.2.2 for implications of attenuation 

wliicli is neglected liere 

3 where /"i o / = 1 and f-^ ^ 1// 




Ionization Parameter [erg cm s 



Figure 1. Fit to temperature- ionization parameter calculated by 
Owen et al. (2010). 



One may also write this at fixed R/Rg to obtain the mass- 
scaling of the density, namely n oc M^^ . Thus, provided we 
know the form of /(C), the density at the sonic surface is also 
fixed a prion through the temperature-ionization parameter 
relation, which provides a density-sound-speed coupling. 



2.1.2 Consequences of Density-Sound-Speed Coupling 

The density-sound-speed coupling found at the sonic surface 
in optically thin heated X-ray heated flows provides some 
very powerful constraints on the flow properties of X-ray 
photoevaporative winds, and implies that the structure and 
temperature of the underlying disc has little effect on the 
resulting flow. 

This arises since in any steady state flow mass-flux must 
be conserved; therefore, in this case since we are able to 
described the mass-flux at the sonic surface a priori with 
reference to the X-ray physics alone, we have determined 
the mass-flux in the entire flow without reference to the disc 
properties. 

This picture is very different to pure EUV (HoUenbach 
et al. 1994) or FUV (Gorti & HoUenbach 2008,2009) flows 
where the underlying density structure in the launch region 
is crucial to the mass- loading of the flow. Unlike EUV/FUV 
flows where the flow will adjust to changes in the underlying 
disc's structure, an X-ray flow will adjust the underlying disc 
structure to 'feed' the flow at the required rate, provided the 
sonic surface remains optically thin to the X-rays responsible 
for the heating. 



2.2 Estimates of the total mass-loss rate 

Armed with the knowledge that it is the X-ray physics, and 
not the disc structure that sets the overall mass-loss rates 
we are free to integrate Equation [3] over the sonic surface to 
estimate the total mass-loss rates. We can use Equations|4]&: 
[7]to substitute for and Cs in the integrand. Furthermore, 
if we change the variable of integration in Equation [3] from 
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cylindrical radius R to escape temperature we find: 

/ 8i-DTr2 r°° 71-1/2 

M^{M.,Lx)= LxJ^^ / dr,,, (8) 

« 8 X lO-'-* (^-jT^ A Mo yr-i (9) 

\ 1 X 10^" erg s^ / 

(where the second result comes from evaluating the in- 
tegral out to an escape temperature of lOOOK, or a ra- 
dius of C'(100)AU for a solar mass star). This mass-loss 
rate is in good agreement with the numerical result of 
6.7 X 1O~^M0 yr~^ obtained from the hydrodynamic sim- 
ulations of Owen et al. (2011b) and is also consistent with 
the approximately linear scaling between mass-loss rate and 
X-ray luminosity found in the simulations. Furthermore, the 
results indicates there is no explicit dependence between 
stellar mass and mass-loss rate (although the M* — Lx scal- 
ing introduces a strong positive implicit scaling). We also 
derive this scalings in more detail through the existence of 
similarity solutions presentecQ in Appendix A. 



2.3 Importance of Numerical Results 

In Sections 4 and 5 we present hydrodynamical simulations 
that have been designed to test and refine the ideas set out 
above. However, the role of the numerical simulations is not 
merely one of validation: although Equation is of great 
utility in predicting the integrated mass-loss from the sys- 
tem (and can be used to derive the mass loss rate per unit 
area at the sonic surface), it does not - without knowledge 
of the detailed streamline topology below the sonic surface 
- allow one to calculate the mass loss rate per unit area 
from the disc itself (Su,(i?)). For example the form of the 
integrand in Equation [5] suggests that the mass- loss is dom- 
inated in the lower temperature regime (2000-5000K), but 
these stream-lines connect the sonic surface to the disc sur- 
face in a non-uniform way, such that there is not a one 
to one correspondence between the mass-flux profile at the 
sonic surface and the mass-flux proflle from the disc. This 
results in a mass-loss profile that peaks at smaller radii in 
the disc, rather than the large radius suggested by the cal- 
culations at the sonic surface. Although this may appear 
to be a detail, it turns out that, when photoevaporation is 
combined with secular disc evolution, the pattern of disc 
clearing is rather sensitive to the form of Su,(_R). In par- 
ticular, it is the form of (R) that determines where the 
gap first opens in the disc, the extent to which the system 
undergoes a period of 'photoevaporation-starved accretion' 
prior to gap opening and hence the surface density profile 
of the disc outside the gap. This latter then determines the 
rate at which the disc is evaporated once the inner disc is 
completely cleared out. This is of obvious importance if one 
wants to compare the predictions of photoevaporation the- 
ory with the observed incidence of transition discs with holes 
of different sizes (see Section 8.1). 



* Note that the scalings presented in Appendix A do not rely 
on the approximate expression for the temperature at the sonic 
surface contained in EquatioiHI 
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Figure 2. Ionization structure of the X-EUV wind around a 
O.7M0Star with an X-ray luminosity of 2 X lO'^^erg s~^, the 
streamlines are equispaced in 2% of the cumulative mass-loss rate 
and the magenta contour shows the region of 98% ionization frac- 
tion (i.e. roughly the penetration depth of EUV photons). This 
shows that an X-ray driven wind is itself entirely optically thick 
to EUV photons and prevents them from reaching the disc surface 



3 THEORY OF UV & X-RAY 
PHOTOEVAPORATION 

We now turn to considering the situation where X-ray heat- 
ing is combined with ultraviolet heating, considering in turn 
the case the ionising (EUV) and non-ionising (FUV) con- 
tinua. 



3.1 EUV and X-ray heated flows 

The radiative transfer calculations of Ercolano & Owen 
(2010) demonstrated why, in the presence of an X-ray flow, 
EUV heating can be neglected. Figure[2]shows the ionization 
structure of an X-ray wind subject to combined X-ray and 
EUV irradiation, with the magenta contour indicating the 
point at which the flow develops a significant neutral fraction 
and beyond which the EUV cannot penetrate. This is easily 
understand if we consider the Stromgren radius: 

/ >1>. \ ^''^ / n N -2/3 

where n = 10* cm^'^ is the density at the base of the flow for 
an X-ray luminosity of Lx = 2 x 10^" erg s~^. Now, Equa- 
tion [7] (or Appendix A) implies the density scales linearly 
with X-ray luminosity (n oc Lx)- If we furthermore assume 
that the EUV and X-ray luminosity are both coronally pro- 
duced (e.g. Alexander et al. 2005), so the luminosities scale 

— 1/3 

linearly we find Vs oc Lj^ . This means that even at the 
lowest X-ray luminosities (~ 10^* erg s~^) the EUV pene- 
tration is still too small to penetrate the X-ray flow. Even 
if we adopt the strongest EUV and the weakest X-ray lu- 
minosities quoted in the literature (i.e. Lx = 10^* erg s~^, 
$, ^'^""^ s-^ Giidel et al. 2007; Alexander et al. 2005) we 
find that the EUV fiux cannot reach the surface of the disc 
from which the bulk of the X-ray wind is launched. There- 
fore, EUV irradiation can be neglected when X-ray radiation 
is also present. 
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(a) X-ray dominated photocvaporation 
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(b) FUV dominated photocvaporation 

Figure 3. Figure showing a cartoon sketch of photocvaporation 
from a disc exposed to UV and X-ray irradiation. The flow begins 
in the FUV heated region, but at some height the flow transitions 
to the X-ray heated region and then passes through the sonic 
surface in the X-ray heated region. 



3.2 FUV and X-ray heated flows 

Gorti & HoUenbach (2008,2009) considered the combined ef- 
fect of UV and X-ray heating in setting the hydrostatic disc 
structure and then used their static structures to predict 
resulting flow rates. In the case of FUV heating, the under- 
lying heated disc structure determines the properties at the 
sonic surface and a variety of solutions, whether static or 
flowing, asymptote to the same (nearly hydrostatic) struc- 
ture at small z. Although the results of Gorti & HoUenbach 
still represent the state of the art in terms of the thermal and 
chemical structure of X-ray and FUV irradiated discs, the 
implications of these studies for the resulting mass loss rate 
are currently unclear. In Section 5 we explore this further 
by using a simple parametrisation of FUV heating in our 
radiation hydrodynamical simulations (although we caution 
this should not be regarded as a proper treatment of the 
FUV heated flow). Here we discuss the general theoretical 
framework for combined FUV and X-ray heated flows. 

X-ray heating is effective up to a total hydrogen col- 
umn of around 10^^ cm~^. For models using depleted dust 
in the disc atmosphere (as in the models of D'Alessio et 
al. 1998,1999,2001), FUV heating extends to significantly 
higher columns. This means that travelling along a line from 
the star to a point in the disc one first traverses gas that is 



exposed to both X-rayfl and FUV, then to FUV only and 
then to neither. Within ~ 100 AU the X-rays dominate the 
heating in the regions that are exposed to both X-rays and 
FUV. Therefore, in this region the sequence (as one proceeds 
away from the star) is X-ray heated, FUV heated and then 
heated by neither - in reality heated by coupling to the dust 
- (e.g. Gorti & HoUenbach, 2009). 

Though we expect this sequence of X-ray heating fol- 
lowed by FUV heating, there are still two topologically 
distinct possibilities which we illustrat^fl in Figure 3. Fig- 
ure |3(a)| is a schematic illustration of the case where the 
fiow in the FUV heated region never exceeds its local sound 
speed for any of the streamlines shown. The sonic transition 
occurs in the X-ray heated region. In this case we can use 
the previous arguments to find that the sound speed at the 
X-ray heated sonic surface is simply given by Equation [3] 
to order unity as before: in other words the temperature at 
the sonic surface is simply a function of cylindrical radius. 
Since this also fixes the density at the sonic surface this im- 
plies that the mass fiux is just a function of R and X-ray 
luminosity. In other words, whatever the detailed physics of 
the heating and cooling processes in the X-ray dark region, 
they must combine to drive a flow that is able to deliver the 
required mass flux to the X-ray sonic surface. We therefore 



expect that if the geometry is as shown in Figure 3(a) the 



photocvaporation rate is entirely controlled by X-ray heat- 
ing. 



On the other hand. Figure 3(b) shows three streamlines 
where the two at larger radius undergo a sonic transition 
within the FUV heated region. In this case the flow has to 
self-adjust to deliver the required conditions at the FUV 
sonic surface - in other words it is the structure of the FUV 
heated region that controls the mass flow rate. Although in 
principle gas flowing out in this FUV heated wind could be- 
come X-ray heated at some point, this region would not be 
causally connected to the flow base, the two regions being 
separated by a shock. However, it is unlikely that an FUV 
driven wind would ever be significantly heated by X-rays 
from the central star. This is because we require the condi- 
tion at the sonic surface in the FUV heated region to satisfy 
Equation m Given that FUV heating does not generally at- 
tain temperatures in excess of ~ 1000 — 2000K (Gorti & 
HoUenbach 2004; Bruderer et al. 2009), this implies that the 



schematic FUV sonic surface shown in Figure 3(b) must lie 



at radii of > 100 AU where X-rays are incapable of heating 
the gas to greater temperatures than the FUV. 

We therefore conclude that FUV driven winds are only 
likely to be significant at radii > 100 AU. If the flow topol- 
ogy at smaller radii is as shown in Figure |3(a)| then the 
photocvaporation rate is set entirely by the X-rays. The 
only situation in which FUV heating could control mass 
loss from much smaller radii is if the streamlines crossing 
the FUV sonic surface at > 100 AU actually originated 



By X-rays we are referring to the X-ray photons that provide 
significant thermal heating with a maximum energy of l-2KeV 
and not those at higher energies that can penetrate larger columns 
but provide insignificant heating (see Ercolano et al. 2009b). 
^ Note the similarity of this discussion to the analysis of fiow 
from proplyds in the case of combined irradiation by FUV and 
EUV radiation fields: see Johnstone et al 1998 
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at small radii Naturally, we need to perform a radiation- 
hydrodynamic simulation with FUV heating to assess this 
possibility (Section 6). 



4 NUMERICAL MODELS OF 

PHOTOEVAPORATING DISCS 

As discussed in Section 12.31 numerical radiation- 
hydrodynamic simulations are required to obtain accurate 
mass-loss rates and importantly the mass-loss profile. In 
this section we extend the parameter space investigations 
by Owen et al. (2010,20ffb) to cover the entire range 
of X-ray luminosity, mass & inner hole size expected for 
discs around low-mass stars, and thus test and calibrate 
the theory discussed above. Owen et al. (2011b) already 
investigated the explicit dependence of mass-loss rate on 
X-ray luminosity and inner- hole size around a O.7M0 star. 
Therefore, in this work we are only concerned with extend- 
ing the models to lower mass and varying disc structures, 
and we have performed simulations of photoevaporating 
discs around 0.1 & 0.7 M© star. 

4.1 Numerical Method 

The numerical method is similar to that described in Owen 
et al. (2010), and is briefly described here for completeness. 
The simulations were performed using a modified version of 
the ZEUS-MP/2 cod^ll (Elayes et al. 2006), where the modifi- 
cations are described in detail in Owen et al. (2010). X-ray 
heating is included using a temperature-ionization param- 
eter relation up-to a critical column of A'^ = lO^^cm"^; at 
columns larger than this value the gas temperature is fixed 
to the dust temperature. We start from the hydrostatic disc 
calculations of D'alessio et al. (2001) that only include dust 
heating; this disc structure is then irradiated with X-rays 
which drive a photoevaporative wind that evolves to steady 
state. We also use spherical co-ordinates in order to match 
the symmetry of the solution at large radius. 

4-1.1 Primordial Discs 

In order to test the mass-scaling, along with the predicted 
X-ray luminosity scaling in different mass regimes, we per- 
form 5 simulations of photoevaporation for a primordial disc 
around a O.IM© star, evenly spaced in logarithmic X-ray 
luminosity from logj^p Lx = 28.3 to 30.3, which covers the 
range of observed X-ray luminosities around M-type stars. 
We use a grid with 288 logarithmically spaced cells in the 
radial direction and 144 cell equally spaced in the angular 
direction. The grid extends from r = 0.03AU (~ 0.03-Rg) to 
r = 46AU (> 457?g) and thus covers the range where we 
expect photoevaporation to be important. 

In order to consider the effect of varying the disc struc- 
ture we run three further simulations of the original calcu- 
lation performed in Owen et al. (2010), i.e. a disc around 
a O.7M0Star with an X-ray luminosity of 2 x 10"^" erg s~^: 

^ Note: In order to increase computational efficiency we have 
moved from the scalar ZEUS-2D code to the MPI ZEUS-MP/2, which 
parallelises efficiently on up-to 64 processors. 
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Figure 4. Mass-loss rate plotted against X-ray luminosity, the 
open squares show the simulations results for a O.7M0Star ob- 
tained in Owen et al. (2011b), the filled circles show the simula- 
tion results for a O.lMQStar described in this work. The solid line 
shows the fit to the mass-loss rates used for a population-synthesis 
study conducted in Owen et al. (2011b), the dashed line shows 
the order of magnitude estimate obtained in Section 3. 

firstly we increase and decrease the dust temperature by a 
factor of -\/2, since the dust in the region that interacts with 
the X-ray heated region is optically thin. These changes in 
temperature effectively correspond to a change in the stel- 
lar luminosity by a factor of 4. We also consider a simula- 
tion with a disc mass of 10% of the original calculation i.e. 
Md = 2.6 X 10"^Mq. 

4-1.2 Discs with inner holes 

Since the low stellar mass simulations presented here have 
lower disc masses than the higher mass counterparts de- 
scribed in Owen et al (2011b), we have to take into account 
of the fact that the density structure of the underlying disc 
may be being eroded during the hydrodynamical simula- 
tion on a time-scale less than the time required to set up 
a steady state wind (this latter being ~ the sound crossing 
time to the edge of the grid; see discussion in Alexander et 
al 2006a). Such erosion is not necessarily 'realistic' since we 
do not model the viscous effects in the disc which may re- 
supply material to the disc inner edge. In order to explore 
this effect we consider three classes of simulation: one with a 
small inner hole (0.7 AU) where the disc is not significantly 
depleted over the duration of the simulation, one starting 
with inner hole radius of 3 AU that was allowed to erode 
over the duration of the simulation and one with radius of 
5 AU where the inner edge is 'held' at 5 AU by artificial 
re-supply (i.e by re-setting the X-ray dark 'base density' to 
its original value every time-step) to obtain a steady-state 
simulation. For each of these cases we model two X-ray lu- 
minosities (logLx = 29.3,30.3). 

4.2 Simulation results: primordial discs 

In Figure |4] we show how the mass-loss rate scales with both 
mass and X-ray luminosity; the results confirm the theoreti- 
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Figure 5. Figure showing the flow topology around a 0.1 (left panels) and 0.7 (right panels) MQStar for X-ray luminosity of 2 X lO'^^erg 
s~^ and 2 x 10^" erg respectively; the results for the O.TMQStar are taken from Owen et al. (2011b). The top panels show the density 
structure with streamlines at 5% intervals of the cumulative mass-loss rate. The lower panels show the temperature structure of the flow, 
where the solid line represents the sonic surface. We note (in agreement with the theoretical model presented in Section 2) that the 
mass-loss profile and temperature are independent of stellar mass when considered in terms of a radius scaled by rg, but the density 
structure scales as MT^ . 



cal predictions described in Section 2.2 to first order, namely 
no explicit mass-dependence on the photoevaporation rate 
and furthermore a mass-loss rate that scales approximately 
linearly with X-ray luminosity. 



In Figure [5] we compare the flow morphology for an X- 
ray luminosity of 2 X 10^^ erg s~^ around the 0.1 (left pan- 
els) and O.7M0 (right panels) stars. The possible existence 
of self-similarity described in Appendix [X] implies that the 
temperature of the flow should be fixed by the value of the 
effective potential, thus to first order the temperature in the 
fiow should be fixed by the mass of the central star and thus 
be scale free when considered in terms of a radius scaled 
by Tg . The lower panels of Figure [5] show that to first order 
this the case. Furthermore the top panels show the density 
structure in the fiow and the streamlines at 5% intervals of 
the cumulative mass loss rate. Despite small differences the 
similarity of the right and left hand panels demonstrate that 
the simulations are close to being self-similar. Moreover, the 
density is lower at the higher mass (at given R/Rg) by a 
factor consistent with the scaling predicted in Section 

2 and Appendix A. We thus find that the simulations are in 
good agreement with the self-similarity arguments presented 
in Appendix A. 



4-2.1 Dependence on disc structure 

As first indicated by the numerical models calculated in 
Owen et al. (2010) with adiabatic discs, the disc structure 
has very little bearing on the final mass-loss rate. In Sec- 
tion 2 we developed an argument to explain this insensitivity 
to underlying disc structure, invoking the necessity to sat- 
isfy particular conditions on the X-ray sonic surface. Here 
we test this idea further by exploring the winds produced 
by more realistic variations in disc structure. 

These correspond to variations in the stellar luminos- 
ity and underlying disc mass. Given the theoretical descrip- 
tion of the fiow predicts that the sonic surface controls the 
mass-flux, in Figure [6] we show the temperature at the sonic 
surface as a function of cylindrical radius. The black line 
shows the escape temperature as a function of radius. We 
see that beyond ~ Rg the models are close to the escape 
temperature, in line with the expectation of Equation 4. 
At smaller radii, the escape temperature rises more steeply 
with decreasing radius than the models. We can simply un- 
derstand this change in behaviour by noting that it occurs 
at a temperature of around 10*K, which is where the ionisa- 
tion parameter temperature relation flattens (see Figure 1). 
At higher temperatures the temperature is rather insensitive 
to ionisation parameter. Thus, whereas at larger radius, the 
ionisation parameter adjusts to produce the temperature at 
the sonic surface demanded by Equation 4, at small ra- 
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Figure 6. The temperature at the sonic surface as a function of 
cylindrical radii (scaled by Rg) for 3 models with two different 
stellar masses (0.1 & 0.7Mq) and X-ray luminosities of (2 X 10^^ 
& 2 X 10^*^ erg s~^). The open squares show the original O.TMq, 
2 X lO'^'^erg s~^ calculated by Owen et al. (2010), the black crosses 
are those of the O.1M02 X lO^^erg s~^ and the open circles are for 
the 0.7Mq , 2 X 10^*^ erg s~^ model calculated with at disc depleted 
in mass by a factor of 10 from the original calculation. The black 
line shows the temperature expected from the order of magnitude 
calculations performed in Section 2. We note all primordial disc 
calculations have sonic temperatures that follow similar loci. 

dius the temperature is more or less fixed. In this case, the 
Parker wind sonic transition corresponds to an almost con- 
stant spherical radius. This transition can be seen in the 
topology of the sonic surface at R ~ Rg in Figures 5 & 8. 
We further demonstrate this behaviour in Figure 7 which 
plots, for a variety of models, the effective potential on the 
sonic surface as a function of cylindrical radius - evidently 
all models evolve so as to satisfy a common relation on the 
sonic surface, regardless of the mass and temperature of the 
underlying disc. 

Figure 8 illustrates the flow topology for the models 
with various disc temperatures and disc mass. Evidently 
these are qualitatively similar but differ in detail. For ex- 
ample, the sonic surface is at lower height and the density 
at intermediate height above the disc plane is reduced when 
the disc mass is lowered. This however has a minimal ef- 
fect on the disc mass loss profiles: even a change in height 
of the sonic surface by 20 — 30% has a small effect on the 
effective potential on the sonic surface and hence hardly af- 
fects the associated temperature, density and mass flux. 
Figure [7] shows that for most radii the value of the effective 
potential is as expected. The dip at i? < 5AU represents 
the region described above that is unable to reach the re- 
quired escape temperature and thus the sonic surface oc- 
curs at a constant height and the z ^ R assumption breaks 
down. This can be understood in terms of the model de- 
scribed in Section 2: the sound speed at the sonic surface 
scales as ct ~ ^teff ~ GMt/2R[2 - (z/R)"^], so that pro- 
vided z ^ R, order unity variations in height correspond 
to small variations in the effective potential (see Figure [7|). 
Thus the sound speed (and hence density and mass flux) 
at the sonic surface are insensitive to such changes in the 
height of the sonic surface. Since the mass-flux is prede- 
termined then the sub-sonic flow structure adjusts itself to 
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Cylindrical Radius [AU] 

Figure 7. The value of the effective potential {h'^ /2R? — GA/*/r 
- note it is shown as a positive value for simplicity) at the sonic 
surface plotted against cylindrical radius for various underlying 
disc structures. The blue points are for the original calculation of 
Owen et al. 2010, the black points represents the mass depleted 
disc described above, the red points show the disc with increased 
dust temperature while the green points show the disc with de- 
creased dust temperature as described above. The dashed line 
shows the prediction from the model described in Section 2. 

feed the sonic surface at this required rate. This results in 
a flow topology that is fairly independent of the temper- 
ature or density in the disc, even though it is ultimately 
these quantities that determine the actual location of the 
sonic surface. As expected we flnd that the total mass-loss 
rates and proflles are fairly invariant to dramatic changes 
in disc structure: the total mass-loss rates for the modifled 
discs range between 0.9 — 1.7 x 1O~*M0 yr~^compared to 
the original 1.4 x lO'^M© yr"^(Owen et al. 2010) for the 
same X-ray luminosity. 

4-2.2 X-ray attenuation 

The basic theoretical outline described in Section 2 assumed 
that the heating at the sonic surface is dominated by opti- 
cally thin X-ray photons. In the numerical models we also 
adopt an optically thin X-ray heating model with a hard 
column cut-off at A'^ = 10^^ cm~^. In order to assess the con- 
sequences of these simplifications, we have performed addi- 
tional simulations where the ionization parameter is instead 
given by: 



In the case of a monochromatic X-ray source, this would 
simply replace the X-ray flux with its value attenuated by 
the appropriate optical depth at that energy. However, in 
reality we have a broad range of energies, and it is less clear 
how to attenuate the flux - in practice the attenuation factor 
will reflect that for the photons that dominate the heating 
locally. Use of the ionisation parameter-temperature rela- 
tion sidesteps the need to follow the propagation of photons 
of different energy and we thus do not have the information 
to model the attenuation properly. However, we can explore 
the consequences of attenuation by toy models in which the 
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Figure 8. Density and velocity structure for photoevaporative winds with same stellar parameters (M, = 0.7Mq, Lx = 2 x 10^'^ erg 
s~^) but varying disc structures. From left to right; Original calculation by Owen et al. (2010); a reduction in dust temperature by \/2; 
a reduction in the disc mass by a factor of 10. All figures shows very similar velocity and density structures in the wind. 
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Figure 9. Temperature structure of the wind from a star with 
an X-ray luminosity of Lx = 2 X 10^'^ erg s~^ and a stellar 
mass 0.7Mq. The solid magenta line shows the sonic surface, dot- 
dashed line a column of 10^^ cm~^, the dashed line a column of 
10^^ cm~^, the dotted line a column of 10^" cm~^ and the solid 
black line a column of 10^^ cm~^. Cold regions at columns < lO'^^ 
cm~^ indicate regions where the dust temperature is higher than 
the temperature given by X-ray heating. 



optical depth of unity is fixed to three characteristic columns 
N = 10^°, lO^S 10^^ cm"^ The models with X-rays attenu- 
ated at columns densities of 10^^ and 10^^ cm~^ show little 
difTerence in structure and mass-loss rate to that performed 
by Owen et al. (2010). However, the model with attenua- 
tion at a column density of 10^° cm~^ is unable to heat the 
outer regions of the disc and the flow is restricted to the in- 
ner -R < 10 AU, resulting in a significantly reduced mass-loss 
rate. We can understand this result by inspecting Figure 9, 
which superposes contours of equal column density to the 
source on a temperature map in the case of a model without 
attenuation. We see that the contour with a column to the 
star of 10^^ cm~'^ lies close to the sonic surface and thus 
(given the arguments that we have given here about the im- 
portance of the sonic surface in setting the mass flux) we 
expect that attenuation at column densities of 10^^ cm~^ or 
above will not greatly affect the mass flow. 



We now need to ask whether the extreme attenuation 
model (i.e. that with an attenuation column of 10^" cm~^, 
which does suppress the mass-loss) is physically reasonable. 
Such an attenuation column would be appropriate to pho- 
tons of energy ~ 0.1 keV. Thus to argue that lO'^" cm~^ 
is the relevant column density is to argue that photons of 
higher energy (i.e. > 0.1 keV) are minor contributors to 
heating gas to the temperatures (3000 — 5000K) in the re- 
gion of the sonic surface that dominates the mass-loss. 

The most convincing argument that the higher energy 
photons (> 0.3 — 0.4 keV) are dominating the heating to 
temperatures < 5000K is provided by the exercise conducted 
in Owen et al. (2010). In this exercise, the wind structure 
of a high density {Lx = 2 x 10^° erg s~^) unattenuated 
steady state wind (obtained from a hydrodynamic calcu- 
lation using an ionization parameter approach) was fed in 
to the MOCASSIN Monte Carlo radiative transfer code (Er- 
colano et al. 2003,2005,2008), which self-consistently follows 
the attenuation and propagation of X-ray photon packets as 
a function of energy. The results of this exercise, shown in 
Figure 9 of Owen et al. (2010), confirmed that the agreement 
between the parametrised temperatures used in the hydro- 
dynamics and the temperatures generated by MOCASSIN is 
excellent. We infer from this that it is higher energy photons 
(> 0.3 — 0.4 keV) that are mainly responsible for heating the 
region of the sonic surface that dominates the mass-flux to 
the required temperatures (< 5000K). 

4.3 Simulation Results: inner hole discs 

In Figure [10] we show the mass-loss rates as a function of 
inner-hole size for the low mass stars (O.IMq) compared 
with previous results of the higher mass-stars. 

As in the simulations reported in Owen et al. (2011b) 
we flnd that the mass-loss rate is approximately indepen- 
dent of the inner hole size. Figure [TO] also shows that the 
steady-state simulations (squares) give identical results to 
the eroding simulations, giving confldence that the mass- 
loss rates obtained from the eroding simulations are accurate 
over the range of radii considered. The simulations show as 
expected that the mass-loss rate scaling is also linear with X- 
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Figure 10. Figure showing the mass-loss rate as a function of 
inner hole size for discs around a O.lMQStar. The filled points are 
those for an X-ray luminosity of 2 X 10"^'' erg s~^ and the open 
points are those for an X-ray luminosity of 2 X 10^^ erg s~^. The 
crosses are the mass-loss rates calculated from the eroded inner 
hole models, while the circles are the moving average calculated 
from these eroding mass-loss rates and the squares are the mass- 
loss rates from simulations which are maintained in steady-state. 
The dashed line shows the mass-loss rates for primordial discs 
with the same X-ray luminosity. 



ray luminosity in agreement with the theoretical arguments 
presented in Section 2 and with previous simulations. 



4-3.1 Understanding the Inner hole Scalings 

Examination of the streamline topology for various inner 
hole radii (over a factor of ~ 2) suggests that the topology 
of the innermost streamline is approximately scale free (i.e 
with an over all radial scale that varies with 7?hoic); in conse- 
quence the variation of temperature and velocity along this 
streamline (as a function of distance along the streamline 
scaled to -Rhoic) is approximately the same in each simula- 
tion (see Figure [11}. This encourages us to construct a scal- 
ing argument even though in the Appendix we show models 
with different -Rhoie cannot be exactly self-similar. There is 
no a priori reason one would expect the temperature struc- 
ture of these flows to be identical: this only follows if we use 
the numerical result that the streamline topology is fixed. 
We argue in the Appendix that the fact that the relation- 
ship between the ionization parameter and temperature is 
not a simple power law implies that the conditions along two 
topologically identical streamlines can only be scaled to each 
other if the variation of temperature with (scaled) distance 
along the streamline is the same in both cases. Then the 
mass-flux scales as n{l)A{l). For fixed temperature Equa- 
tion [7] implies n oc R^^i^ and the stream-bundle area scales 
as A oc flhoie- Therefore, along this streamline the mass-flux 
is fixed and does not scale with inner hole radius. 

The origin of a scalable innermost streamline is likely to 
originate from the fact that unlike primordial discs (where 
the thermal pressure is dominant in determining the stream- 




Figure 11. Figure showing the streamline topology and den- 
sity, sound speed and velocity along the sub-sonic section of the 
streamline. The streamline topology is shown in the top left hand 
panel scaled in terms of -Rholci with the sonic surface shown as the 
dotted line. Then bottom left panel show the density X R^^^^, the 
top right show the sound speed and the bottom right show the 
velocity as a function of distance along the streamline normalised 
to ij^ojg. The black lines are for a hole of 14 AU, the blue line 
for a hole of 21AU and the blue for an inner hole of 30 AU. 



line), in this case the centrifugal force is important since 
material that is stripped from the rim of an inner hole ini- 
tially moves inwards. As the effective potential is scalable in 
terms of iihoie then it follows that (in the limit of negligible 
pressure) the streamline topology would also be scalable in 
terms of -Rhoic. Provided the total mass-flux is dominated 
from the inner hole region - as the simulations indicate - 
then the total mass-loss should scale in a similar fashion to 
the mass-flux on the innermost stream-bundle, resulting in 
a total mass-loss rate that is approximately independent of 
the inner hole radius. 



5 NUMERICAL TESTS: INCLUDING 
'FUV-TYPE' HEATING 

In Section 3 we argued that at radii < 100 AU the flow mor- 



phology is likely to resemble that shown in Figure 3(a) where 
the mass loss rate is set entirely by the conditions at the X- 
ray sonic surfaces and where any FUV heating deeper must 
produce a flow that self-adjusts to supply the correct mass 
flux at the X-ray sonic surface. In order to test whether the 
flow topology resembles that shown in Figure |3(a)| we per- 
form some simulations that roughly aim to include the effect 
of extra-heating below the X-ray dominated region. If FUV 
heating was important to the calculation of photoevapora- 
tion rates then one would require a detailed calculation of 
X-ray and FUV heating and cooling, along the lines of com- 
bining the X-ray algorithm with an FUV algorithm (e.g. 
Richling & Yorke, 2000). While such a calculation would be 
necessary to predict line emission from the warm atomic and 
molecular layer, we are here interested in the simpler ques- 
tion of whether a plausible modification of the structure at 
the base of the flow due to the effect of FUV heating will in- 
fluence the flow at the X-ray sonic surface. Therefore, a full 
X-ray -|-UV calculation is unnecessary and beyond the scope 
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of this work and we instead employ a simple prescription for 
FUV heating at the base of the flow. 

We calculate the FUV flux throughout the disc/wind 
system in an identical manner to Richling & Yorke (2000), 
including opacity only from dust extinction (e.g. we ignore 
opacity arising from the ionization of Carbon and Sulphur, 
thus overestimating the flux and gas temperature). We then 
turn the attenuated FUV flux into a temperature using a 
parametrisation of the previous FUV calculations of Rich- 
ling & Yorke (2000) and FUV heating calculation of Gorti 
& HoUenbach (2008): 

/ J. \ 0.35 / X -0. 

Tpuv = 1400K ( A ( ^ A 

V300 erg s-i cm-2y V5 x lO^cm-V 

(12) 

where /fuv is the attenuated FUV flux. This parametri- 
sation reproduces the temperature structure of the region 
directly irradiated by FUV in Richling & Yorke (2000) 
(their model A) and the temperature structure listed for 
the 'launching point' of the calculations of Gorti & HoUen- 
bach (2008). We perform three simulation with a large FUV 
luminosity of 1 x 10''^ erg s^^ where we vary the FUV ex- 
tinction, and choose a depleted (from ISM type dust) value 
of the FUV extinction of 300 cm^ g~^ used by Richling & 
Yorke (2000) and a extinction that is depleted by a further 
factor of 10 to 30 cm^ g~^. Furthermore the calculations of 
Gorti et al. (2009) indicate that more depleted discs drive 
more efficient flows (as in their calculations the mass-loss 
remains constant as the accretion rate - and thus FUV lu- 
minosity - falls). Therefore, we perform a calculation using 
the low mass disc described earlier {Md = 2.6 x 10~^Mo) 
and use an FUV extinction of 2 cm^ g~^ indicative of a very 
heavily settled disc atmosphere (D'Alessio et al. 1998,1999). 
At every point in the disc/wind system we choose the heat- 
ing mechanism (X-rays, FUV or dust) that produces the 
highest temperature, where we note the dust temperatures 
used already include the effects of FUV heating from accre- 
tion luminosity as well as the chromospheric contribution 
(D'alessio et al. 2001). Furthermore, as we attenuate the 
FUV flux there is no need to apply a column cut-off as in 
the X-ray case. Thus provided the FUV gas temperature 
is higher than the X-ray or dust temperature at any point 
in the grid (even those with Nh > 10^^ cm~^) it can be 
heated by the FUV radiation. The FUV extinction is fixed 
throughout the calculation domain: this is reasonably valid 
as Owen et al. (2011a) showed that photoevaporative winds 
are capable of entraining small dust particles, which will 
dominate the opacity at FUV wavelengths. We emphasise 
that this is no attempt to calculate the self-consistent ther- 
mal structure of a disc heated by FUV/X-rays, but purely an 
investigation into the influence of FUV type heating at the 
base of an X-ray heated flow. We demonstrate this further 
by running calculations in which we vary the parametrisa- 
tion; while this of course has no physical basis it is purely a 
simple way to demonstrate that are our results are not be- 
ing driven by our detailed assumptions about FUV heating. 
We thus calculate an extra 4 models for the maximum FUV 
penetration model, where we vary the power law indices (in 
Equation 12 in all four possible combinations) of the flux 
and density variation by factors two in both directions. 

Finally, we investigate whether the X-ray to FUV lu- 
minosity ratio has any bearing on the results (particularly 



any transition from the flow topology shown in Figure 3(a) 
to Figure 3(b) when the X-rays are weak compared to the 
FUV luminosity). We thus fix the FUV luminosity at the 
large value of 1 x 10^^ erg s~^ and systematically lower the 
X-ray luminosity to the lowest values observed ~ 10^* erg 
s-i. 



5.1 Results 

Our calculations using this extra 'FUV-type' heating does 
not affect the mass-loss rates obtained (for Lx = 2 x 
W^°, Lfuv = 1 X 10^^ erg s"^): we find a value of 

1.2 X 1O~^M0 yr~^for the initial calculation with the same 
depleted FUV extinction coefficient as Richling & Yorke 
(2000), a value of 1.1 x 10~**Mq yr~^for the calculations 
with the very depleted FUV extinction coefficient, and a 
value of 1.3 x 10"**Mq yr"^for the mass depleted disc. In 
all simulations we find as expected that the sonic surface 
occurs in the X-ray heated regions and we therefore confirm 
the simple ideas presented in Section 2, as well as demon- 
strating that the photoevaporation rates are controlled by 
the X-ray physics alone. 

In Figure[T2]we show the density and velocity structure 
of the maximum FUV penetration model (i.e. the mass de- 
pleted disc). This shows as expected that the sonic surface 
occurs in the X-ray heating region, while the FUV heated 
flow is confined to a layer between the bound dust heated 
disc and the X-ray heated flow (i.e. between the two black 
lines in Figure I12|) . We further find that even though the 
density and temperature structure of the 'FUV heated re- 
gion varies between model runs, in all cases the thermal 
structure of the sonic surface is as shown in Figure [S] There- 
fore, the FUV layer must be adjusting its structure in order 
to 'feed' the X-ray heated region with the correct mass-flux, 
and satisfy the hydrodynamical conditions at the sonic sur- 
face. Therefore, extra heating below the X-ray heated region 
has little effect on the photoevaporation rates in the range 
7? ~ 1 - lOOAU. 

The structure of the system shown above does allow in 
principle for some additional mass loss due to FUV heating 
if some of the gas between the two black lines flows radi- 
ally outwards through the FUV heated wedge rather than 
feeding material to the X-ray sonic surface. If such a 'wedge 
flow' (see schematic depiction in Figure [13)1' did exist then 
(given the temperature of the FUV heated gas) it would not 
undergo a sonic transition until a radius of order 100 AU, 
i.e. close to the outer boundary of the grid. We find that gas 
does leave the grid sub-sonically in the FUV heated region. 
However, due to the use of outflow boundary conditions on 
this grid boundary, it is impossible to tell if this flow is real 
or an artefact of outflow boundary conditions (which are 
only exact for super-sonic flow). Therefore, we are unable to 
further hypothesise on the possible existence of such FUV 
'wedge' flows. 



5.1.1 Effect of X-ray luminosity 

In Figure [14] we show how the mass-loss rates vary with 
decreasing X-ray luminosity, compared to the mass-loss rates 
expected from pure X-ray models. For Lx > erg s"'^ 
the flow resembles that of Figure 12: the X-ray sonic surface 
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Figure 12. Velocity and Density structures of winds including an 
extra 'FUV like heating mechanism, the arrows show the velocity 
in the wind and the black lines show regions where the extra 
heating mechanism is dominating over the X-rays and dust and 
the magenta contour shows the sonic surface. The model shown is 
the depleted disc with the extremely low FUV extinction allowing 
for maximum penetration of the underlying disc. 




Dust Heated Disc 



R~100 AU 

Figure 13. Cartoon of an X-ray flow with an FUV wedge flow 
underneath. 



extends out to cylindrical radii > 50 AU and sets the mass 
loss rate from the disc. At lower X-ray luminosities (for a 
fixed Lfuv of 10^^ erg s~^) the subsonic FUV region puffs 
up vertically and drives the X-ray sonic surface and upper 
region of the FUV heated region almost vertical at a radius 
of ~ 10 AU. This means that the predominant mass loss 
in the disc becomes at this point dominated by sub-sonic 
flow from the FUV region. As noted above, we do not trust 
the mass loss rates that we obtain in the case that material 
leaves the grid sub-sonically so that the mass loss rates at 
the low Lx end of Figure 14 should not be taken as a good 
measure of the total mass-loss rates. Thus considerably more 
investigation is required to describe photoevaporation in the 
limit of low X-ray to FUV luminosity ratios. 

Instead, the chief result of Figure 14 is that the flow 
morphology and mass-loss rates are X-ray dominated pro- 
vided that Lx/Lyw sJ 0.01. While stars with X-ray lu- 
minosities < 10^^ erg s~^ only represent a small fraction 
^ 5 — 10% of the population of young stars (Giiedel et al. 
2007; Owen et al. 2011b), the role of combined FUV/X-ray 
heated flows will certainly be important for these objects if 
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Figure 14. Mass-loss rates arising in the X-ray heated flow from 
models that include the FUV heating parametrisation, plotted 
as a function of X-ray luminosity. The dotted line indicated the 
mass-loss rate expected from a pure X-ray model. See text for a 
discussion of the reliablity of the results in the limit of low X-ray 
luminosities. 



the FUV luminosity continues to remain large (~ lOOLjc) 
throughout the disc's lifetime. 



6 FINAL CLEARING OF THE DISC 

In Figure [T5] we show the evolution of a simulation with 
Af. = 0.1Mq& Lx = 2 X 10^" erg s"^ in which the disc 
hole is eroded during the course of the simulation (and for 
which the mass-loss rate during the evolution is shown by the 
crosses in Figure 10) . During the first 300 years of evolution, 
the photoevaporative flow slowly strips off gas from the inner 
edge of the disc. During this stage, there is a radially thin 
layer of bound X-ray heated gas {Tx ~ 200 — 400K) on the 
inside rim of the disc and the flow direction of the escaping 
gas is radially inwards. The origin of this layer is fairly easy 
to understand: as the fixed streamline topology suggested by 
the simulations fixes the density and structure in the fiow, 
the flow cannot by construction consist of a flxed column 
density. It turns out that the column density of the flow 
is much less than 10^^ cm~^. Therefore, the X-ray heated 
region must contain not only the flow but also a warm bound 
X-ray heated region, that is separated from the flow by a 
contact discontinuity. As the hole grows, the radial width of 
the layer of warm bound X-ray heated gas grows in size in 
order to still maintain the critical column density for X-ray 
heating of lO^'^ cm"'^. Once the hole size grows to i?hoic ~ 
lOAU, the radial width of this region becomes comparable 
with the vertical scale height of the X-ray heated gas. At 
this point, the topology of the X-ray heated gas just within 
the rim changes from nearly radial inflow to one where there 
is a strong (roughly sonic) expansion velocity normal to the 
disc plane. This is seen in Figure [15] as the plume at the 
inner edge of the outer disc which, in the lower panels, rolls 
over to envelope the disc. The importance of this change of 
flow topology is that it allows gas to be evacuated out of 
the line of sight between the X-ray source and the disc rim. 
This thus increases the penetration of X-rays in the mid- 
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O.lMQStar. The colour map shows the density structure and the arrows indicate the velocity structure. Each panel indicates a snapshot 
of the flow structure. 



plane. The panels at 353.6 years and 392. years illustrate 
the formation of a low density bubble of heated gas between 
the expanding plume and the residual disc. One can also 
discern that this bubble is expanding radially into the disc 
and compressing the disc's inner rim. We term this effect 
'thermal sweeping' and note that it operates on a roughly 
dynamical time-scale (i.e. ~100 years at lOAU) which is not 
much more than the free expansion time of X-ray heated 
gas. We emphasise that this evolution is much faster than 
what occurs in the early stages of clearing inner hole discs, 
when the erosion rate is limited by the fact that the flow is 
radial and hence the evaporating gas continues to provide a 
degree of shielding to the disc's inner rim. 

This 'thermal sweeping' was not seen in the simula- 
tions calculated in Owen et al. (2010,2011b) for discs around 
a Q.TMq star. We can understand this difference in that 
our simulations have a constant initial disc to star mass ra- 
tio and the previous simulations therefore had higher disc 
masses and higher mid-plane densities. This implies that 
the critical column for X-ray absorption was attained over a 
correspondingly smaller radial length scale; thus, at a given 
radius, the ratio of radial thickness to vertical scale height 
of the X-ray heated gas is smaller in the simulations with 
higher disc mass. These simulations therefore remained in 
the regime where the stripped gas moved radially inwards; 
we expect that (if we had the computational resources to 
calculate inner hole models around higher mass stars with 
very large radii) we would at some stage have entered the 
regime where the heated layer became radially thick and 
that 'thermal sweeping' would ensue in this case also. 

When this X-ray warm, bound region is small it must 
be in dynamical balance with the photoevaporative flow and 
the dust heated disc. The time-scale on which it can obtain 
this equilibrium is the sound crossing time of the region 
tac ~ A/cs, where A is the width of this region. When 
A/H <^ 1 (where H is the vertical scale height), this re- 



gion can adjust to dynamical equilibrium quickly and the 
disc/flow system is stable, even though the pressure in the 
region is significant. However, at the point when A ~ 
this region cannot rapidly obtain dynamical equilibrium, 
and therefore will expand vertically on the same time-scale. 
This vertical expansion can clearly be seen in the time se- 
ries plots, and as it expands the mid-plane density drops, 
resulting in higher temperatures. This in turn accelerates 
the vertical dispersal and a runaway ensues. Eventually, a 
high pressure region develops at the edge of the dust heated 
disc that pushes the remaining disc material away rapidly 
on the dynamical time-scale. 

In order to model this process we use the fact that 
the warm bound X-ray heated region must be in dynam- 
ical equilibrium with both the photoevaporative fiow and 
the dust heated disc. The pressure is largest where bound 
X-ray heated region meets the dust heated disc, this region 
is highly sub-sonic and the requirement of dynamical equi- 
librium is simply one of pressure equilibrium. Thus one may 
write: 



ksnxTx = Pdust 



(13) 



where nx is the density of the bound X-ray heated region. 
Now we can use Equation ll3l to find A: 



Nx 
nx 



(14) 



where Nx is the X-ray penetration depth (i.e. Nx = 10^^ 
cm~^). Thus the requirement for this process to begin A = 
H, simply becomes: 



^dust(^) ^ 



NxkgTx 
H{R) 



Y.{R) < V2TTfimhNx 



Tx 

Tdust 



(15) 



(16) 
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Figure 16. Plot of the LHS (solid- pressure of dust heated disc) 
and RHS (dashed-critical pressure of X-ray heated region) of 
Equation 1151 to estimate the inner hole radii at which thermal 
sweeping begins. The stars show the radii the simulations indicate 
this process begins. Black represents a disc around a O.TMq star 
while blue shows a disc around a O.IMq (shown in Figure [TSj. 



Substituting standard values into this equation we find: 

(17) 

Thus this process will begin when the mid-plane pressure 
of the dust heated disc has dropped to a sufficiently small 
value, mainly through a reduction in the mid-plane density. 
This can be achieved either on account of a low disc mass or 
because the hole has grown to sufficiently large radius. This 
is shown in Figure [161 where in order to test the applicability 
of Equation 15 we also perform a re-simulation of model 
F from Owen et al. (2010) (a disc with an initial hole of 
20 AU, an X-ray luminosity of 2 x 10^" erg s~^ and mass 
of O.7M0), but the mid-plane pressure has been reduced 
to the point where Equation [15] predicts this process will 
occur when the disc hole size modestly exceeds its initial 
value. This re-simulation shows similar morphology to that 
shown in Figure [15] and the disc is completely cleared via 
the 'thermal sweeping' mechanism described above as the 
hole reaches ~ 25 AU (shown in black on Figure I16p . The 
remaining disc mass just before this process sets in is around 
5-10 Jupiter masses, and therefore a small fraction of the 
disc's initial mass is cleared in this process: roughly ~ 80% 
of the disc's initial material in accreted onto the star, ~ 
10% lost in a photoevaporative wind and the final ~ 10% 
disperses through this process of 'thermal sweeping'. 

However, it is worth mentioning that these simulations 
were designed to study photoevaporative flows in the opti- 
cally thin X-ray heated gas. On the other hand, once 'ther- 
mal sweeping' is initiated, the flow becomes sonic close to 
the location where X-ray heating is initiated and in this case 
a simple column cut-off will not accurately reproduce the 
thermal structure of this region. Furthermore, as the gas 
is considerably cooler (hundreds to thousands rather than 
thousands to tens of thousands K as in the photoevapora- 



tion calculations) we cannot neglect the effect of molecular 
cooling, along with the added effect of FUV heating. This 
latter depends on the dust properties at the mid-plane of 
the disc at large radii: although the X-ray heated region is 
probably opaque to FUV radiation in the case of dust that is 
ISM-like, the FUV may penetrate beyond the X-ray heated 
region in the case of dust that is significantly depleted. The 
necessary calculations (including X-ray and FUV radiation 
along with dust evolution) are certainly beyond the scope 
of this work. We here restrict ourselves to stating that our 
simplified calculations suggest that the disc enters a 'ther- 
mal sweeping' phase (which clears the residual outer disc 
on roughly a dynamical time-scale) once the disc enters a 
region when the mid-plane pressure is sufficiently small. 

Finally, while the instability sets in on dynamical time- 
scales, simple energetic calculations suggest it may take ~ 
10 dynamical time-scales (at the inner edge) to unbind the 
outer regions of the disc at the lowest X-ray luminosities. 
This is somewhat longer than our simulations would predict 
in this case, since the simulations - by imposing radiative 
equilibrium - neglect the effect of adiabatic cooling. Hence, 
adiabatic cooling may be important during the final stages of 
'thermal sweeping', but not during the initial stages (onset 
of the instability) as the amount of material contained in the 
bound X-ray heated region is small compared to the total 
disc mass. However, 10 dynamical time-scales (at the inner 
edge, ~ 10* years) is still considerably faster than any other 
gas removal process operating at this stage. 



7 IMPLICATIONS FOR DISC EVOLUTION & 
DISPERSAL 

The ultimate goal of the photoevaporation theory is to be 
able to predict disc evolution and lifetime based on prop- 
erties of the underlying systems. At its most basic level 
we have shown throughout this work that it is only stellar 
X-ray luminosity that sets the photoevaporation rate and 
hence is the major factor in determining the disc lifetime. 
We have argued that the vigorous photoevaporation rates 
(10~* — lO"'^^ Mq yr~^) are an inescapable consequence 
of X-ray irradiation (provided X-rays reach the disc's sur- 
face, an issue discussed in detail in Section I8.2|l , Further- 
more, we showed that variations in disc structure and the 
inclusion of EUV and FUV radiation fields have little effect 
on the fiow structure and derived mass-loss rates provided 
Lx/Lfvv > 0.01. 

The major consequence of this model is of course that 
discs around stars (of the same mass) with higher X-ray lu- 
minosities lose their discs first. An obvious implication of 
such a result is that disc-less stars will have on average a 
higher X-ray luminosity than disc-bearing stars at the same 
age. This has been observationally recognized for a while 
(e.g. Preibisch et al. 2005) and Owen et al. (2011b) showed 
that the X-ray photoevaporation model is quantitatively 
consistent with these observations. Ingleby et al. (2011) has 
attempted one of the first studies to investigate the role of 
high energy irradiation on disc dispersal; unfortunately the 
available X-ray data for the study was incomplete in both 
X-ray luminosity and mass, meaning it was impossible to 
disentangle the effects of stellar mass and X-ray luminosity 
on disc dispersal and the result was inconclusive. Now we 
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have fully developed the photoevaporation model down to 
lower masses, we would have been hopeful that it would be 
possible to make predictions of disc evolution as a function 
of stellar mass. However, as described in Section 18.11 the 
disc's lifetime also depends on how the viscous time varies 
with stellar parameters and speculation about its variation 
with stellar mass is well beyond the aims of this paper. It is 
worth noting though that as discussed in Section 8.1, sen- 
sible assumptions lead to the conclusion that disc lifetime 
is not a strongly dependant function of stellar mass and is 
likely to decrease only weakly with increasing stellar mass, in 
agreement with the observations (e.g. Ercolano et al. 2011b). 

T.l Transition Discs 

It was the observations of the frequency of 'transition' discs 
(Strom et al. 1989; Skrutskie et al. 1990; Kenyon & Hart- 
man, 1995) that lead to the development of the photoevapo- 
ration mechanism. Since the original observations, the num- 
ber of detected discs that are consistent with the original 
characteristics of a transition disc (i.e. a large drop in opac- 
ity at short wavelengths, while returning to primordial levels 
at longer wavelengths) has increased. However, there is no 
clear observational or theoretical consensus as to what a 
'transition' disc actually is and it is becoming clear that the 
population of observed 'transition' discs are not a homoge- 
neous sample (e.g. Alexander & Armitage, 2009; Owen et al. 
2011b). The observations now span a range that includes ac- 
creting, non-accreting objects and those with small amounts 
of gas and dust inside the opacity 'hole'; it is unlikely that 
one mechanism is responsible for all objects (Cieza et al. 
2010). Some 'transition' disc samples are now being split 
into 'transition' and 'pre-transition' populations (e.g. Es- 
paillat et al. 2010; Furlan et al. 2011). 

The three main mechanisms that can account for the 
formation and structure of this class of objects are: photoe- 
vaporation, planet formation and grain growth. Owen et al. 
(2011b) showed that a large fraction (> 50%) of solar type 
'transition' discs were consistent with the X-ray photoevap- 
oration scenario. In the X-ray photoevaporation scenario a 
'transition' disc goes through two stages of evolution. Once 
the gap opens, the dust in the residual inner disc is rapidly 
(i.e on a time-scale much less than the viscous evolution time 
of the gas) depleted as a result of inward migration due to 
gas drag. The loss of dust (and its associated opacity) from 
the inner disc would result in the object being classified as a 
transition disc even while it still possesses an inner gas disc 
with its associated accretion. This phase lasts approximately 
~ 10% of the observable disc's lifetime. Once the inner disc 
has entirely drained onto the central star, the inner edge of 
the outer disc is directly exposed to X-ray irradiation and 
the outer disc is then eroded to large radius. During this 
phase a centrifugal barrier prevents material accreting from 
the inner-edge onto the star and the disc would appear as a 
non-accreting 'transition' disc with a large inner hole. The 
structural properties of observed 'transition'/'pre-transition' 
disc are currently more detailed than the photoevaporation 
model can predict. This must await the development of cou- 
pled dust and gas evolution models that are correct when 
the dust-to-gas ratio is no longer small. Therefore, at this 
stage we do not distinguish between 'transition' and 'pre- 
transition' discs. Although discussed in Section 7.2 the in- 
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Figure 17. The inner hole radius and accretion rate of low-mass 
objects classified as 'transition' discs by varies authors: cyan tri- 
angles (Brown et al. 2009); blue filled triangles & squares (Cieza 
et al. 2010); black open squares (Ercolano et al. 2009); black filled 
squares (Espaillat et al. 2008,2010); black open circles (Hughes 
et al. 2009, 2010); red filled circles (Kim et al. 2009); magenta 
filled diamonds (Merin et al. 2010); green filled star (Najita et al. 
2010). Where error-bars or ranges are listed they are shown. Any 
object classified as non-accreting it is placed on the plot with an 
accretion rate of 10~^^Mq yr~^. The grey region represents the 
region of parameter space where the objects are consistent with 
gaps created by photoevaporation. See text for more details. The 
blue and black dashed regions represent those regions for a 0.3 
and IMq star respectfully. 



elusion of 'thermal sweeping' significantly changes the ac- 
creting to non-accreting fraction predicted by previous pho- 
toevaporation models, to a population that will be a major- 
ity of accreting transition discs. 

The observational data however also contains examples 
of transition discs with large inner holes (i?hoic ii 20AU) and 
large accretion rates (M ~ 1O~^M0 yr~^) that are simply 
impossible to create through any kind of photoevaporation 
model. It has been suggested that these objects could have 
their large accretion rates set by FUV ionization and MRI 
activation at the inner edge of the hole (Perez-Becker & Chi- 
ang, 2011); in this case the inner hole (or gap) is created by 
the tidal effect of a planet within the hole and dust filtra- 
tion at the inner edge of the outer disc is invoked to reduce 
the opacity of the inner disc. However, such models are at 
an early stage of development and is not clear what kind of 
planetary system would be required to match the accretion 
rates and spectral energy distributions in these systems. 

Using the photoevaporation theory developed in the 
previous sections we can delineate the region of transition 
disc parameter space within which we expect objects created 
by photoevaporation to lie for stars with A/, < 1.5MQ(This 
exercise is similar to that undertaken by Owen et al. 2011b 
except that we have now extended the predictions to differ- 
ent stellar mass ranges and added new observational data 
points). Since the radial scale of the photoevaporative fiow 
scales linearly with mass, the radius at which the gap opens 
will obviously scale linearly with mass too. Furthermore, 
since the gap opens when the accretion rate falls below the 
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photoevaporation rate, then stellar mass effects are not im- 
portant in the determining the accretion rates of transition 
discs other than through the implicit stellar mass depen- 
dence of the X-ray luminosity (and hence photoevaporation 
rate) . 

In Figure [T71 we show the region in which photoevapo- 
rating transition discs are expected, over-plotted with obser- 
vationflaf discs classified as transition discs: cyan triangles 
(Brown et al. 2009); blue filled triangles & squares (Cieza et 
al. 2010); black open squares (Ercolano et al. 2009); black 
filled squares (Espaillat et al. 2008,2010); black open circles 
(Hughes et al. 2009, 2010); red filled circles (kirn et al. 2009); 
magenta filled diamonds (Merin et al. 2010); green filled star 
(Najita et al. 2010). We note that as Cieza et al. (2010) do 
not fit for an inner hole radius but rather list the Spitzer 
band of the inner hole. Therefore, we split the sample into 
two: the objects with a Aturnoff ^5 5.8 fim, the inner hole is 
conservatively set to be less than 5AU, for Aturnoff = 8 
we again conservatively set an inner hole of lOAU (since 
the inner-hole must occur in the temperature range A = 8 
- 24 /im), we add that these are safe overestimates based 
on simple temperature structures of protoplanetary discs. 
It is also worth noting that at very low masses (~ O.IM©) 
Ercolano et al. (2009) showed that current observations can- 
not determine the presence of holes smaller that ~ 1 AU, 
due to a lack of contrast between the disc and the stellar 
photosphere. 

The photoevaporation region is constructed by using 
the fact the photoevaporation profile is self-similar when 
the radius is scaled in terms of Rg. (i.e. if the gap opens 
at lAU around a IMQStar it would open at O.IAU around 
a O.lM0Star) and we use a X-ray luminosity scaling of 
Lx oc M^^"^ (Preibisch et al. 2005, Guedel et al. 2007) to 
extend the region calculated by Owen et al. (2011b) to the 
full mass range. Again this shows that the observations can 
be separated into two groups, those with small inner-holes 
< 20AU and small accretion rates < 10^* M0 yr^'^that are 
consistent with a photoevaporative origin, and those discs 
with large accretion rates and large holes that cannot have 
a photoevaporative origin. 

7.1.1 Role of stellar mass 

Kim et al. (2009) noted that 'transition' discs may show cor- 
relations between disc properties and stellar mass. Within 
the X-ray photoevaporation framework, such correlation will 
be primarily driven by the variation of the X-ray luminosity 
with stellar mass {Lx oc m5^; Preibisch et al. 2005, Guedel 
et al. 2007). In order to assess the role stellar mass may play 
in driving correlations between 'transition' disc properties, 
we have extended the population synthesis study performed 
in Owen et al. (2011b) to lower (O.IM©) and higher (IM0) 
masses, where we adopt an initial disc mass that scales lin- 
early with stellar mass. Furthermore, as the scaling of the 
initial viscous time with stellar mass is unconstrained (see 

* Several other types of astrophysical objects are found to enter 
these samples (e.g. AGB star or debris discs): in these cases we 
follow the definitions of the authors and only plot objects that 
are classified as protoplanetary discs with a gap/hole in its dust 
disc. 
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Figure 18. Panels showing the predicted regions occupied by 
transition discs created by photoevaporation. The shaded region 
represents the parameter space occupied by accreting transition 
discs (i.e. those with M* > 1 X IO^^-'Mq yr~^) , the dashed line 
represents the maximum radius a transition disc may reach (even 
when a non-accreting transition disc) before 'thermal sweeping' 
sets in. 

Section 8.1), we preform runs where t^, oc Af» and t^, oc Af^^, 
although the results shown are only mildly sensitive to the 
chosen scaling. 

In Figure [18] we show the predicted regions that the 
inner hole radius and accretion rates of transition discs oc- 
cupy as a function of X-ray luminosity and stellar mass, both 
for accreting transition discs (A/, > 1 x 10~'^^Mq yr"'^), and 
non-accreting transition discs before 'thermal sweeping' sets 
in. In all cases they show strong positive correlations, driven 
by the radial scaling of the photoevaporating region with 
mass {Rg oc Af*), and the positive correlation between mass 
and X-ray luminosity. The shaded regions are calculated by 
randomly sampling in time the disc during the transition 
phase, and represent the area populated by all transition 
discs through the entire population's evolution. Thus, the 
correlations may vary with time for individual clusters, as 
only a certain fraction of discs will be in the transition disc 
phase at that point in time. 

7.2 Relic Discs 

As discussed in Section [6] we have uncovered a further phase 
in disc dispersal, i.e. the rapid removal of disc material due 
to what we term 'thermal sweeping'. This phase begins at 
the point that the bound X-ray heated gas just inward of 
the inner disc rim becomes geometrically thick and when 
this gas starts to fiow predominantly vertically, thus allowing 
further X-ray penetration in the disc mid-plane. The popula- 
tion synthesis calculation carried out in Owen et al. (2011b) 
of solar type stars (Af» — 0.7Mq) did not include this effect 
and therefore predicted the existence of long lived relic discs, 
for a small fraction of the total population. These were discs 
with very large inner holes R > lOOAU and low photoevap- 
oration rates; for objects with the lowest X-ray luminosities 
these relics were estimated to remain for lOMyr and Owen 
et al. (2011b) argued that such discs might well correspond 
to some sources that had previously been classified as young 
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Figure 19. Cumulative distribution of inner hole radii (solid line) 
and disc mass (dashed line) at which the thermal sweeping process 
starts for the disc models calculated in the Owen et al. (2011b) 
population synthesis model (O.7M0Stars). 

(gas free) debris discs. This new thermal sweeping efTect now 
provides a mechanism for the removal of these relic discs on 
a rapid (dynamical) time-scale (~ 10^~^ years). 

In order to investigate at what point we expect the ther- 
mal sweeping to set in for the case of discs around higher 
mass stars (the calculations reported in Section [6] were for 
stars of mass Q.IMq) we apply the pressure criterion derived 
in Equation 1151 since this corresponds to the point at which 
the radial width of the X-ray heated bound gas is ~ H. 
We follow the population synthesis approach of Owen et al. 
(2011b), i.e. we assume that all discs evolve according to a 
viscous similarity solution with fixed parameters and where 
the only quantity that varies from star to star is the X-ray 
luminosity which sets the photoevaporation rate. Then each 
star-disc system is evolved under the combined effects of vis- 
cous evolution and photoevaporation, following each system 
through the phase of gap opening and outer disc clearing. 
We identify the point that 'thermal sweeping' sets in with 
point that the maximum mid-plane pressure in the viscous 
disc falls below the value given by Equation 15. The varia- 
tion of this radius with X-ray luminosity is rather modest, 
given that X-ray luminosity affects the radial profile of the 
clearing disc only through the mild effect of the degree of 
depletion of the disc during the phase of photoevaporation 
starved accretion (Drake et al. 2009; Owen et al. 2011b). 
Consequently our population synthesis model predicts that 
in the higher mass stars, thermal sweeping should set in at a 
narrow range of radii around 40 AU and disc masses < lOMj 
(see Figure 19). 

Furthermore, the identification of this gas expulsion 
method significantly alters the ratios of accreting/non- 
accreting 'transition' discs expected, since discs spend less 
time in the non-accreting phase than previously expected. 
In Figure [2171 we plot (as a function of X-ray luminosity) the 
fraction of the disc's lifetime that the disc spends as a tran- 
sition disc (crosses) and as a non-accreting transition disc 
(dots). This indicates that discs with high (> 10^" erg s~^) 
and low (< 10^^ erg s'^) X -ray luminosities spend the ma- 
jority of their transition phase as accreting objects. Discs 
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Figure 20. Fraction of the disc's total lifetime that the disc 
presents as a transition disc (crosses) and a non-accreting transi- 
tion disc (dots) as a function of X-ray luminosity. 

with intermediate X-ray luminosities spend similar amounts 
of time in the accreting and non-accreting phases. Thus over 
the entire population of discs, one would expect most transi- 
tion discs to appear as accreting. We can understand the rel- 
ative lack of non-accreting transition discs at both high and 
low X-ray luminosity as follows. At low X-ray luminosity the 
disc accretes most of its mass onto the star over its longer 
lifetime, so that the mid-plane pressure is low once it en- 
ters the accreting phase and thermal sweeping sets in rather 
promptly. At high X-ray luminosity the process of photo- 
evaporation starved accretion (Drake et al. 2009; Owen et 
al. 2011b) is vigorous and once the disc enters the clearing 
phase the outer disc is very depleted of gas (compared to the 
result of pure viscous evolution). Again the mid-plane pres- 
sure is low and the disc enters the thermal sweeping phase 
quickly. 

The reason that such thermal sweeping was not seen in 
our previous hydrodynamic simulations of these more mas- 
sive stars (Owen et al 2010) is that our hydrodynamical 
calculations of discs with inner holes used input disc profiles 
that were simply those of primordial discs truncated at a 
given value and did not take into account the fact that the 
outer disc profile would have been modified by viscous evo- 
lution and photoevaporation. The discs modelled thus had 
higher pressures than discs that had been self-consistently 
evolved and were thus less prone to thermal sweeping. This is 
in contrast to the sequence of disc clearing for a 0.1 M© star 
(see Figure 12) in which the evolution of the disc profile is 
followed self-consistently in the hydrodynamic calculation. 
We emphasise that our estimate that the disc is rapidly 
cleared in higher mass (0.7 M©) stars once the hole grows 
to 40 AU is very approximate given that it is not based on 
self-consistent hydrodynamical simulation of disc clearing in 
this case. 



8 DISCUSSION 

In this work we have indicated that one of the inescapable 
consequences of X-ray heating of a protoplanetary disc is 
a vigorous photoevaporative flow. The mass-loss rate is ef- 
fectively determined by the X-ray luminosity only (when 
Lx/Lfuy > 0.01), with disc structure and stellar mass only 
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having weak effect on the resultant flow structure. This is 
easily understood in terms of the model presented in Sec- 
tion 2, where the conditions at the sonic surface are deter- 
mined by the X-ray physics alone, so the disc below has 
to simply adjust to feed the sonic surface with the correct 
mass-flux. The consequence of such large mass-loss rates 
are that photoevaporation will be competing with accretion 
over much of the disc bearing life of young stars. Stars with 
higher X-ray luminosities will begin losing their discs first, 
followed by the lower X-ray luminosities stars at late times. 
Such an evolution is consistent with the observation that 
disc-less stars are systematically more X-ray luminous than 
stars still surrounded by discs (Neuhauser et al 1995; Stelzer 
& Neuhauser 2001; Flaccomio et al. 2003; Preibisch et al. 
2005). 

We have shown it is extremely unlikely that the EUV 
irradiation field plays any role in setting the photoevapo- 
ration rates: in any scenario where the EUV can reach the 
disc then so can the X-rays and the resulting X-ray driven 
wind is optically thick to the EUV. We have also set out 
arguments (backed up by simple toy FUV heating models) 
that FUV heating becomes important only below the X- 
ray sonic surfacqj- Only when the X-ray luminosity drops 
to sufficiently low (Lx/ipuv <0.01) values does the FUV 
begin to shut off the X-ray flow, and may possibly drive a 
flow itself. This result raises the intriguing possibility that 
what determines the evolution of a disc (and presumably its 
ability to form planets) is the star's intrinsic X-ray luminos- 
ity, rather than any variations in initial conditions when the 
star/disc system was created (Armitage 2011). 



8.1 Disc lifetimes as a function of stellar mass 

Since we have shown that the photoevaporation rate is in- 
dependent of stellar mass (at fixed Lx) then, given some 
knowledge of the statistical dependence of X-ray luminosity 
on stellar mass and of the scaling of disc viscous parameters 
with stellar mass, we can examine the predicted relationship 
between disc lifetime and stellar mass. To this end we adopt 
the 'null' disc evolution method (Ercolano & Clarke 2010; 
Owen et al. 2011b), i.e. we assume that discs evolve in a pure 
viscous manner until the accretion rate equals the photoe- 
vaporative wind rate. At this point the disc is cleared rapidly 
and we define this point as the disc lifetime. For a disc in 
which v oc R the viscous similarity solutions (Lynden-Bell 
& Pringle, 1974) teU us: 

M4t) = M40)h+—] (18) 

where ti, is the disc's initial viscous time-scale (which, 
for the assumed linear scaling of viscosity with radius) is 
itself just proportional to the disc initial radius (-Ri). Thus 
we can simply write the disc life-time (rd), when t > ti, as: 

Ta(xM,{Of^''M-''/h, (19) 



^ This is not to say that FUV heating is then irrelevant to setting 
the structure of the warm disc at intermediate heights (Gorti & 
HoUenbach 2008,2009), but merely that is not the determinant of 
the mass-loss rate. 



Now using the fact A/d(0)/A/* (0) oc and Mw oc Lx we 
find 

Td OC MdiOf^^L-^^h'J' (20) 

Observations and theoretical work that suggest the 
discs enter this viscous phase after a self-gravitating phase 
suggest that Md{0) oc A/*, Furthermore X-ray observa- 
tions (Guedel et al. 2004, Preibisch et al. 2007, Albacete 
Colombo et al. 2007) indicate that Lx oc M^^^~^''^ , adopt- 
ing Lx oc Mt^^ we find that: 

Td OC M-'^h'J'' (21) 

Unfortunately the scaling of (or equivalently Ri ) with 
stellar mass is unknown. Some suggested scalings have been 
motivated by a wish to reproduce the claimed quadratic de- 
pendence of accretion rate on stellar mass (although it is 
likely that this empirical correlation - Natta et al. 2006 - 
may be largely driven by a mixture of sensitivity and selec- 
tion effects: see Clarke & Pringle 2006, TiUing et al. 2008). 
One prescription that yields this quadratic dependence is 
that of Alexander & Armitage (2006), who assume that 
most observed discs are at an age that is much less than 
their initial viscous time-scale. Under this circumstance, the 
accretion rate is simply linked to the ratio of initial disc 
mass to initial viscous time-scale; a quadratic scaling of ac- 
cretion rate on stellar mass then demands that (assuming 
Md{0) oc M») the initial viscous time-scale oc 1/M*. If 
this were true, then Equation [21] would imply that disc life- 
time scaled as M^'^^^ . There are however a range of other 
possibilities based on reproducing the claimed accretion rate 
versus stellar mass relation: e.g. if Md{0) oc M« and ti, is in- 
dependent of stellar mass (as proposed by DuUemond et al. 
2006) then Td oc A/y^, while if we adopt A/d(0) oc M, and 
now assume instead that discs are observed at times » 
then we would require tv oc Af« in which case we again ob- 
tain Td OC MI^^ . 

Clearly our estimates are resulting in a range of depen- 
dences none of which has a strong stellar mass dependence. 
Given the form of Equation 1211 we see that a strong de- 
pendence of dispersal time-scale on stellar mass could only 
result if the initial disc viscous time-scale were an extremely 
strong function of stellar mass. This would appear to be un- 
likely, both because there are no obvious theoretical reasons 
to invoke such a strong mass dependence and also because 
such a strong dependence would imply a large stellar mass 
dependence of the relative numbers of stars observed with 
discs that are a substantial fraction of the mass of the central 
star. 



8.2 Model Uncertainties 

It is worth spending some time dwelling on some of the un- 
certainties and possible sources of error in a direct applica- 
tion of the model to real objects. These can be divided into 
two categories: (i) assumptions that directly affect the cal- 
culated mass-loss rates; (ii) how the mass-loss rates affect 
disc evolution. 
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8.2.1 Uncertainties in mass-loss rates 

The analysis presented in Section 2 suggests that the mass- 
flux through the sonic surface is dominated by gas with 
temperatures around 2000 — 5000K. At these temperatures, 
molecular cooling (which we neglected in the X-ray ther- 
mal calculation used here, which just consider line cooling) 
could in principle be important (Glassgold et al. 2004; Mei- 
jerink et al. 2008; Gorti & HoUenbach 2008,2009). In order 
to assess whether this is in fact the case (given that the 
densities at the sonic point are many orders of magnitude 
lower that those encountered in the preceding disc calcula- 
tions) we need to estimate the relative contributions of the 
atomic line cooling at the sonic surface with that provided by 
molecules. From our calculations, the atomic line cooling, is 
typically ~ 10~" erg cm"^ s"^ at 2000K. Given the density 
at the sonic surface (~ lO'' cm"''), then molecular cooling 
would become important only if the CO abundance (Xco) 
exceeded of 0.1 at 2000K (using Equation C7 of Glassgold 
et al. 2004) which is many orders of magnitude higher than 
expected in discs (e.g. Bruderer et al. 2009). Thus we are 
confident that molecular cooling can be neglected at the 
sonic surface in photoevaporation calculations and hence the 
mass-loss rates will be unaffected. 

While in this work and Owen et al. (2010) we went 
to great effort to indicate that the majority of factors that 
vary with time (e.g. disc structure, FUV/EUV luminosity) 
make very little difference to the derived mass-loss rate, we 
made an important assumption, namely the X-ray flux that 
is incident on the disc is an unvarying function of time. 
At first thoughts this is not a bad starting point since the 
time-averaged (X-ray variability due to fiares are too fast 
to change a photoevaporative fiow - Owen et al. 2011b) X- 
ray output of the star has been shown observationally to 
be fairly constant as a function of time (Guedel et al. 2004, 
Ingleby et al. 2011), crucially what is important to the pho- 
toevaporation model is not the X-ray output of the star but 
the 'soft' {hv < IkeV) X-ray luminosity visible to the disc. 

Firstly, since most X-ray studies are concerned with the 
total X-ray luminosity of the objects (normally taken out to 
hv ~ lOkeV) , not much attention has been paid to the evolu- 
tion of just the 'soft' component of the X-ray spectrum with 
time. While varying the amount of 'hard' X-ray irradiation 
will not affect the heating and hence photoevaporation rate, 
it will have an effect on the translation between observed 
X-ray luminosity of an object and the total 'soft' X-ray fiux 
that reaches the disc. 

A second factor is that various sources of absorbing col- 
umn lie between the stellar corona and disc surface, for ex- 
ample: accretion columns, jets and (magnetic) disc winds 
may be present. Obviously if the absorbing column from 
these process approaches lO'^^cm"^ then the X-rays will be 
screened from the disc and the mass-loss rates derived in 
the previous sections will be an overestimation of the actual 
photoevaporation rate. Thus, simply assuming that the disc 
sees an X-ray source with the same luminosity as the intrin- 
sic X-ray luminosity may be an oversimplification, leading 
to an over estimate of the photoevaporation rates. 

Finally, although we are confident about the magnitude 
of the X-ray photoevaporation, we cannot rule out the possi- 
bility that beyond 100 AU it is FUV photoevaporation that 
dominates the mass loss rate. This is because this is the 



distance at which the X-rays and FUV fields heat the gas 
to comparable temperatures and where also such temper- 
atures are comparable with the local escape temperature. 
Gorti & HoUenbach (2009) have suggested that this may 
be an important effect in truncating the disc at the outside 
(cf. a similar role for external FUV radiation in the case of 
the Orion proplyds Johnstone et al. 1998; Richling & Yorke, 
2000; Adams et al 2004). Such truncation would influence 
the viscous evolution of the disc (Clarke 2007) by imparting 
a different boundary condition from the one adopted in the 
viscous calculations as noted by Gorti et al. 2009. 

8.2.2 Uncertainties in disc evolution 

Uncertainties in disc evolution do not affect the photoevap- 
oration rates for a given stellar mass. X-ray luminosity and 
(where relevant) inner hole radius, given the insensitivity of 
these rates to the structure of the underlying disc. However 
it is the interplay between viscous evolution and photoevap- 
oration that sets the sequence of gap opening and inner disc 
clearing. Viscous evolution sets the time-scale on which the 
gap opens and, by setting the mass remaining in the outer 
disc, also sets the time-scale of outer disc clearing. 

Currently, the magnitude and even origin of the viscous 
processes in discs is not fully determined, nor the radial 
range within which various viscous processes are active (see 
Armitage 2011 for an up-to date review). Therefore, in order 
to determine the evolution of discs we adopt a similarity 
law put forward by (Lynden-Bell & Pringle, 1974) wherein 
v (X R since this crude assumption appears to be broadly 
compatible with observations of disc mass distributions and 
accretion rates in T Tauri stars (e.g. Hartmann et al. 1998). 

Although the broad sequence of evolutionary phases 
outlined above is to be expected for any viscosity prescrip- 
tion, we are aware that there are quantitative effects that 
do depend on the magnitude and the radial dependence of 
the viscosity. In particular, we note that (as described in 
Section 17.11) . the inner disc rapidly becomes optically thin 
following gap opening, as it becomes significantly dust de- 
pleted. Such a change will undoubtedly change the thermal 
and ionization structure of the disc. This in turn can be ex- 
pected to affect the viscosity in the disc with consequences, 
for example, for the predictions we have made concerning 
the relative incidence of accreting vs non-accreting transi- 
tion discs with a photoevaporative origin. 

9 CONCLUSIONS 

In this work we have outlined the basic properties of the pho- 
toevaporation of discs around low mass stars by the X-rays 
from the central star. In this regime, the photoevaporation 
rates are set by the properties in the X-ray flow and are fairly 
insensitive to variations in disc structure and extra heating 
sources e.g. FUV radiation. At the most basic level this in- 
dicates factors that determine discs' lifetimes and evolution 
may stem from the stars' intrinsic X-ray luminosities. Our 
main conclusions are summarised below: 

(i) The photoevaporative mass-loss rates are primarily set 
by the X-ray heating alone, as the sonic surface always oc- 
curs within the X-ray heated region. Extra UV heating pro- 
vides little effect on the determined structure of the X-ray 
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wind and its associated sonic surface geometry and mass-loss 
rate if Lx/ipuv > 0.01. 

(ii) We have derived from first principles the scahng of 
photoevaporative mass loss from primordial discs with X-ray 
luminosity and stellar mass and have also derived analytical 
estimates of the total mass loss rates. The mass-loss rates 
scale approximately linearly with X-ray luminosity and have 
no explicit dependence on stellar mass. These scalings and 
absolute values have been confirmed through a large suite 
of numerical simulations. 

(iii) We do not predict a strong dependence of disc life- 
time on stellar mass for stars less massive than IMq, in 
agreement with observations. Higher mass stars have higher 
X-ray luminosities but also have more massive discs to be 
dispersed. The resulting time-scale for disc dispersal would 
be a strong function of stellar mass only in the unlikely case 
that the viscous properties of discs were an extremely strong 
function of stellar mass. 

(iv) We have explained why, in the case of discs with in- 
ner holes, the photoevaporative mass-loss rates scale linearly 
with X-ray luminosity and is independent of inner hole size 
and have calculated numerical models to calibrate the mass- 
loss rates. 

(v) We find that once the mid-plane pressure of a disc 
with an inner hole drops below a critical value (given by 
Equations 15 & 16) the remaining disc is cleared rapidly 
(on a roughly dynamical time-scale). This process (which 
we have termed 'thermal sweeping') sets in once the layer 
of bound X-ray gas on the inside of the disc rim becomes 
thick (i.e. comparable with its vertical scale height) since 
at this stage the heated gas escapes preferentially normal 
to the disc plane and increases the exposure of the disc rim. 
Such rapid dispersal removes the possibility of the long lived 
'relic' discs hypothesised in Owen et al. (2011b). We find 
that thermal sweeping sets in at around 10 — 15 A.U. in 
the case of low mass stars (O.IM©) and we roughly estimate 
that the corresponding radius is < 50 AU for higher mass T 
Tauri stars (~ 0.7Mq). 

(vi) We find that when we take into account the spread 
in masses and X-ray luminosities of T Tauri stars, we can 
produce transition discs that span a wide range of properties 
(hole radii and accretion rate onto the star). In this way we 
can account for a large fraction (> 50%) of the observed 
transition discs. We cannot, however, account for large holes 
with high accretion rates (see Figure [17)1 . 

(vii) Photoevaporation and 'thermal sweeping' ultimately 
destroy the final disc, but it is still accretion onto the central 
star that removes most of the original disc mass (~ 10% 
Mt): in most cases only ~ 10 — 20% of the original disc 
mass lost through in the processes of photoevaporation or 
thermal sweeping. 
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APPENDIX A: THE THEORY OF X-RAY 
PHOTOEVAPORATION: SCALING RELATIONS 

In this appendix we develop the theory of X-ray photoevap- 
oration from discs with a view to understanding how local 
quantities vary as a function of model input parameters. 
We make use of the fact that in thermal equilibrium the 
temperature of optically thin. X-ray heated gas is a roughly 
monotonic function of the ionisation parameter ^ — Lx /nr'^ 
In qualitative terms, the flow starts highly sub-sonically at 
the base of each streamline and is then accelerated by a 
mixture of pressure gradients and effective gravity, gcfi (i.e. 
the combined effects of the gravity of the central star and a 
centrifugal term). It is assumed that viscous effects are neg- 
ligible over the flow; thus the specific angular momentum of 
each fluid element is conserved, being at the Keplerian value 
at launch (cylindrical radius Rh) with value hj^ = GM^Rh- 
Therefore, 



S=ff = -^r + ^^R 



(Al) 



where f , R are unit vectors in the spherical radial and cylin- 
drical radial directions. 

In a steady state we can write the continuity equation 
in the form: 



dlogp _ dlogtt dlogA 



(A2) 



dl dl dl 

(where I is a coordinate measured along the streamline 
and A represents the area of a streamline bundle, so that 
dlogj4/d/ is simply the divergence of the unit vector along 
the streamline, 1). This allows us to write the steady state 
momentum equation along a streamline in the form: 



dlogti dlogA dlogCs gcft ■ 1 



d; 



dl 



dl 



+ ■ 



(A3) 



Therefore a flow will undergo a sonic transition at the point 
where the right hand side of Equation IA3I is equal to zero. 
In the special case of an isothermal flow, we would recover 
the result that in the absence of an external force (geff = 
0) such a transition occurs at the point of transition from 
converging to diverging flow (the de Laval nozzle solution). 
Furthermore, in the case of spherical outflow under point 
mass gravity, the right hand side of Equation I A3I vields the 
sonic condition of a Parker wind (Parker, 1958). 

However, the present case is much more complicated, 
since we do not know the flow geometry a priori, because 
QeS is non- negligible at all radii and because Cs is not con- 
stant. Although it is possible to solve for u{l) (and hence 
p{l) and Cs{l)) for a fixed streamline topology, there is no 
guarantee that such a solution would satisfy the steady state 
momentum equation perpendicular to the streamlines. We 
can use Equation IA3I to write the condition for the sonic 
surface as: 



2 GM R fg 

2R r fA 



(A4) 
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where 



and 



-geff • 1 

GM/r'^ 



fA = 



dlogA _ dlogc^ 
dl dl 



2/r 



(A5) 



(A6) 



The first two terms equate the temperature at the sonic 
surface to the usual expression for a spherical Parker wind; 
fg takes into account of the different component of the ef- 
fective force along the streamline compared with the Parker 
case (both due to the centrifugal term and the non-spherical 
geometry), while /a takes into account the different diver- 
gence of non-spherical flows and also the fact the flow is 
non-isothermal. In thermal equilibrium the flow time-scale 
is long compared to any thermal time-scales and the deriva- 
tive of the sound speed is small compared to the other terms 
(e.g. Adams et al. 2004). Furthermore, as the centrifugal 
force falls off rapidly with radius, /a and fg approach unity 
quickly (Begelman et al. 1983; Adams et al. 2004). Therefore 
provided the sonic surface is some distance from the base of 
the flow, the sound speed will be of order the Parker value. 

At first look these remarks appear to preclude further 
analytic arguments, but in fact we can make useful progress 
by deriving scaling relations. First we consider the case 
where we have obtained a self-consistent steady state flow 
structure (via a hydrodynamical simulation) for a given set 
of input parameters, and we then change some parameter 
of the flow and obtain a new structure. We can then ask 
whether there are circumstances where one case is simply a 
scaled version of the original case. 

If such a scaling does apply, then every streamline in 
the flrst case can be identifled with a topologically identi- 
cal streamline in the second case, which differs only in its 
overall geometrical scale. Such a scaling also implies that 
the variation of flow variables along each such streamline are 
scaled versions of each other, i.e. we write all variables in the 
form u = uou{l) where / is the distance along the streamline 
scaled to the radius at the base of the streamline (Rt) and 
Mo allows for a simple re-scaling of the velocity between the 
two streamlines. In this case we can write Equation I A3I as: 



CaO 



dlogfi 
dl 



dlogA 
d[ 



dloge^ 
d[ 



GM, geff • 1 



where 



gcfl ■ 



i?3 



(A7) 



(A8) 



Now in order to obtain a consistent, self-similar solution we 
require that the above equation should contain only scaled 
variables (i.e. those denoted with a tilde) and should not 
depend on the parameters Ri,, uq and c^o which vary from 
one simulation to another. Examining the left hand side of 
Equation I A7I we then require that uo = Cso (i.e. the Mach 
number of the flow is the same for the two scaled streamlines 
at given I). 

We can furthermore argue that the individual values of 
Mo and Cbo are the same for streamlines that are a scaled ver- 



sion of each other: if Cso were different then the temperature 
at given / would be different between the two streamlines (by 
a constant factor) and the density would also be necessar- 
ily different according to the parametrisation of temperature 
against ionisation parameter. Since this relation is not a sim- 
ple power law scaling then the density variation along the 
streamlines would not be self-similar, in contradiction with 
our assumption. Therefore, for streamlines to map onto one 
another in a self-similar way, we must have the same velocity 
and temperature structure as a function of I. 

Turning now to the right hand side of Equation IA7I we 
see that in order that this does not contain variables that 
differ from case to case, we must have Ri,o oc M, . Since the 
escape temperature is given by Tesc oc M,/r, this implies 
that Tesc is the same at given I for streamlines that map 
onto each other in a self-similar way, thus the importance of 
rg as the natural scaling of the problem becomes clear. 



Al Primordial discs 

We flrst consider the case of 'primordial' discs on the as- 
sumption that these are approximately scale free. Power 
law variation of mid-plane variables is an expectation of the 
Lynden-Bell & Pringle (1974) similarity solutions and ap- 
pear to be consistent with observational estimates of disc 
surface density proflles which suggest that E oc R~^ (An- 
drews et al 2010). These would therefore appear to be good 
candidates for yielding self-similar solutions when the input 
parameters are varied, as the boundary conditions of the 
problem would be close to self-similar. 



A 1.1 Variation of Lx at constant M, 

We have seen above that strict self-similarity requires that 
Rbo oc A/* . Thus if we vary Lx at flxed M, , then _Ri,o is the 
same and each streamline simply maps onto a streamline 
with the same radius and with the same temperature and 
velocity variation along the streamline. Since the temper- 
ature is a function of ionisation parameter, it follows that 
as one varies Lx, the density scaling along each streamline 
must simply vary oc Lx. The mass flux along each streamline 
is the product of stream-bundle area, velocity and density. 
As the former two are independent of Lx and the density 
scales with Lx, then the mass flux is simply oc Lx. 

Therefore, we conclude that primordial discs which dif- 
fer only in their values of incident Lx should be topologically 
identical, with identical velocity and temperature structure, 
but with a mass weighting on each streamline that scales 
with Lx . This result agrees with the results of the hydrody- 
namical simulations of Owen et al (2011b) in which it was 
found that the mass loss rate scales linearly with Lx. 



A1.2 Variation of M^, at constant Lx 

In this case, the fact that Rto oc M» implies that variation of 
stellar mass should simply change the over-all radial scale of 
the flow (i.e. a streamline rooted at 10 A.U. in a simulation 
with M = IMq should simply map on to a streamline rooted 
at 1 A.U. in the case of a simulation with M — O.IMq). As 
before, we require that the two scaled streamlines should 
have the same variation in velocity and temperature as a 
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function of distance along the streamline normalised to the 
radius of the streamline base, i.e. the same variation of ioni- 
sation parameter. Since ^ = Lx/nr^, it follows that a radial 
scaling oc M* results in a scaling with density of oc M~'^ . 
The mass-loss from the streamline is proportional to the 
product of the density (oc M^^), the local velocity (con- 
stant) and the area of the stream-bundle (oc M^): i.e. the 
mass loss rate along each streamline is the same as the mass 
loss rate over its scaled equivalent. Therefore (provided the 
disc is sufficiently radially extended that it encompasses the 
whole region from which significant mass loss occurs in all 
cases), we conclude that, at constant Lx, the mass fiux in 
the wind is independen^f\ of M« . 



A2 Discs with inner holes 

Here the truncated density profile of the disc introduces a 
particular radial scale (-Rhoie) to each underlying disc . How- 
ever, if one fixes -Rhoie for a given mass, then the above 
arguments for Lx are still applicable. For fixed M« and 
Rhoie then as one varies Lx the density must also vary as 
n oc Lx, resulting in a total mass-loss rate that varies lin- 
early with X-ray luminosity. Furthermore, if one varies M», 
then one can simply map one solution with an inner hole ra- 
dius to the equivalent inner hole around another mass using 
7?hoie oc M* . 

Neither of these cases however describe the situation 
that is normally of greatest interest i.e. the variation of pho- 
toevaporation rate with i?hoie in the case of fixed M« and 
Lx, which would allow us to compute the evolution of the 
mass loss rate as a disc is progressively cleared from the in- 
side out. Since -Rhoic introduces a fixed scale into the prob- 
lem no similarity solution exists than includes Lx, Mt and 
fJhoic and it is not possible to analytically estimate such a 
scaling. Therefore, we defer a discussion of this case to Sec- 
tion 5 where our numerical results allow us to construct a 
posteriori an approximate scaling argument. 



A3 Summary 

We have shown that the fiow from primordial discs is ex- 
pected to be self-similar (when Lx and A/* are varied) and 
predict that the mass loss rate should scale linearly with Lx 
and be independent of M* . Furthermore, we have shown that 
photoevaporation rates should scale linearly with Lx also in 
the case of discs with inner holes and that a disc with given 
inner hole radius can also be scaled to a different hole size 
(oc M.) when the stellar mass is varied. The flow topology 
is however not exactly self-similar in the case of discs with 
different hole sizes and fixed Lx and M*. Given the lack 
of exactly scalable solutions in this case, we will return to 
this issue, using additional input from numerical solutions, 
in Section 5. 



It is of course important to note, that even though there 
is no explicit mass dependence on photoevaporation rates from 
discs, there will of course be a strong implicit variation with mass 
through the variation of X-ray luminosity with stellar mass: ob- 
servations suggest that M* oc ^^'^ (e.g. Preibisch ct al. 2005, 
Guedel et al. 2007). 



APPENDIX B: NUMERICAL FITS TO 
MASS-LOSS RATES 

In this appendix, we provide fits to the total mass-loss rates 
and profiles obtained from the numerical calculations per- 
formed by Owen et al. (2010), Owen et al. (2011b) and 
in this work. The mass-loss profile fits were based on the 
results provided by the on-line 'function finder' provided 
at "http: / / zunzun.com", and their functional form does not 
represent the results of analytic calculation and should be 
treated in such a way. Furthermore, these fits were per- 
formed in such a way as to described the cumulative mass- 
loss rates accurately which is important for the global vis- 
cous evolution. As such the fits are worst at small radius; 
should sensitivity at these radii become important more ac- 
curate fits to this region can be provided upon request. The 
surface mass-loss profiles Em given in Equations IB2I & IB5I 
are provided in normalised form and must be scaled so that 
J 27ri?E„di? yields the required total mass- loss rate, given 
in Equations lBT1fclB4l 



Bl Primordial Discs 

The total mass-loss rate variation with X-ray luminosity is 
shown in Figure 31 where the solid line is described by: 

/AT \ -0 068 ^ ^ s 1.14 

M„ = 6.25x10-^ ( ^ 1 [ ^ , ) Mq yr-^ 



1 Mr. 



VlO^Oerg s-i 



(Bl) 

The mass-loss profile, derived from the solid line in the left 
hand panel of Figure 4 in Owen et al. (2011b) is described 
by: 

Eu,(a; > 7) = lo*-"! '"^i"'^'"' 

X lO'^'l loglo{^)^ + El loglo(^)^+/l iogio(^)-i-9i) 

6ailog(a:)^ 5bilog(a:)^ 4ci log(a:)=' 



2;2 1og(10)'' a;2 1og(10)6 a;2iog(io) 
2 



3dilog{xf 2e^\og{x) 

a;2 1og(10)4 x2 1og(10)3 x2 log(10)2 



X exp 



Vloo/ 



(B2) 



where ai = 0.15138, h = -1.2182, ci = 3.4046, di 
-3.5717, ei = -0.32762, /i = 3.6064, gi = -2.4918 and: 



X = 0.85 



V AU I 



(B3) 



where Eii,(x < 0.7) = 0. Logarithms of the form log(a;) are 
taken taken using the natural base. 



B2 Discs with Inner holes 

As described previously we find that the total mass-loss rate 
is approximately independent of inner hole radius, which is 
approximately described by: 



M„ = 4.8x10" 



IQSOerg s-i 



M© yr ' 
(B4) 

The mass-loss profile, derived from the solid line in the right 
hand panel of Figure 4 in Owen et al. (2011b) is described 
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by: 



a2&2exp(62y) ^ C2d2 exp(d2y) _^ e2/2exp(/2y) 



R 



R 



R 



X exp 



(B5) 



where aa = -0.438226, 62 = -0.10658387, C2 = 0.5699464, 
d2 = 0.010732277, 62 = -0.131809597, /2 = -1.32285709 
and: 



y = 0.95 (7? - i?hoi< 



1 Mp 



(B6) 



where Eu,(j/ < 0) = 0. 
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